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(57) ABSTRACT 

Disclosed herein is a system for rapidly resolving position 
with centimeter-level accuracy for a mobile or stationary 
receiver [4]. This is achieved by estimating a set of param- 
eters that are related to the integer cycle ambiguities which 
arise in tracking the carrier phase of satellite downlinks 
[5,6]. In the preferred embodiment, the technique involves a 
navigation receiver [4] simultaneously tracking transmis- 
sions [6] from Low Earth Orbit Satellites (LEOS) [2] 
together with transmissions [5] from GPS navigation satel- 
lites [1]. The rapid change in the line-of-sight vectors from 
the receiver [4] to the LEO signal sources [2], due to the 
orbital motion of the LEOS, enables the resolution with 
integrity of the integer cycle ambiguities of the GPS signals 
[5] as well as parameters related to the integer cycle ambi- 
guity on the LEOS signals [6]. These parameters, once 
identified, enable real-time centimeter-level positioning of 
the receiver [4]. In order to achieve high-precision position 
estimates without the use of specialized electronics such as 
atomic clocks, the technique accounts for instabilities in the 
crystal oscillators driving the satellite transmitters, as well as 
those in the reference [3] and user [4] receivers. In addition, 
the algorithm accommodates as well as to LEOS that receive 
signals from ground-based transmitters, then re-transmit 
frequency-converted signals to the ground. 

19 Claims, 15 Drawing Sheets 
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SYSTEM USING LEO SATELLITES FOR 
CENTIMETER-LEVEL NAVIGATION 

I. CROSS-REFERENCE TO RELATED 

APPLICATIONS 5 

This application is a continuation of U.S. patent applica- 
tion Ser. No. 09/287,523, filed on Apr. 7, 1999, now aban- 
doned which is a continuation of U.S. patent application Ser. 

No. 09/167,520 filed on Oct. 6, 1998, now abandoned, 
which is a continuation of U.S. patent application Ser. No. 10 
09/045,497 filed on Mar. 20, 1998, now abandoned, which 
is incorporated herein by reference, which claimed benefit of 
priority from U.S. Provisional Patent Application No. 
60/041,184, filed Mar. 21, 1997. 

This invention was reduced to practice with support from 15 
NASA under contract number NAS8-39225. The U.S. Gov- 
ernment has certain rights in the invention. 

II. BACKGROUND 

20 

Conventional satellite-based positioning techniques are 
based on the use of special navigation signals transmitted 
from several navigational satellites. In the global positioning 
system (GPS), for example, a constellation of GPS satellites 
transmit LI and L2 carrier signals modulated with C/A and 25 
P code signals. By measuring these code signals, a user 
receiver can determine its position to an accuracy of several 
meters. 

To determine the user position with higher accuracy, a 
differential technique can be used. A reference receiver 30 
having a known position also measures the code signals and 
calculates its position. The reference receiver then calculates 
a differential correction by comparing its known position 
with this calculated position, and transmits this correction to 
the user receiver. Assuming the user receiver is near the 35 
reference station, it can use the differential correction data to 
improve the accuracy of its position estimate down to 
approximately 1 meter. 

Various proposed techniques provide positioning accu- 
racy on the order of 1 cm. In addition to measuring the code 40 
signals from the GPS satellites, these techniques use carrier 
phase measurements of the signals from the GPS naviga- 
tional satellites. Typically, this carrier phase positioning 
technique uses differential carrier phase correction data from 
a reference station in order to improve performance. There 45 
is a significant difficulty inherent to this technique, however. 
When tracking a carrier signal of a navigational satellite 
transmission, one is able to directly measure the phase of the 
signal, but one cannot determine by direct measurement how 
many complete integer cycles have elapsed between the 50 
times of signal emission and reception. The measured carrier 
signal thus has an inherent integer cycle ambiguity which 
must be resolved in order to use the carrier phase measure- 
ments for positioning. Consequently, much research in the 
art of satellite-based positioning has focused on resolving 55 
these cycle ambiguities in carrier phase measurements of 
GPS satellite signals. 

MacDoran and Spitzmesser (U.S. Pat. No. 4,797,677) 
describe a method for deriving pseudoranges of GPS satel- 
lites by successively resolving integers for higher and higher 60 
signal frequencies with measurements independent of the 
integers being resolved. The first measurement resolves the 
number of C/A code cycles using a Doppler range; these 
integers provide for independent measurements to resolve 
the number of P code cycles, and so on for the L2 and LI 65 
carriers. This technique, however, assumes exact correlation 
between satellite and user frequency standards (i.e., the user 
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requires an atomic clock), and provides no means of cor- 
recting for atmospheric distortions. 

A similar technique, called dual-frequency wide-laning, 
involves multiplying and filtering the L2 and LI signals 
from a GPS satellite to form a beat signal of nominal 
wavelength 86 cm, which is longer than either that of the LI 
signal (19 cm) or the L2 signal (24 cm). Integer ambiguities 
are then resolved on this longer wavelength signal. Since the 
L2 component is broadcast with encryption modulation, 
however, this technique requires methods of cross- 
correlation, squaring, or partially resolving the encryption. 
These techniques are difficult to implement and how low 
integrity. 

Hatch (U.S. Pat. No. 4,963,889) describes a technique for 
resolving integer ambiguities using measurements from 
redundant GPS satellites. Initial carrier-phase data is col- 
lected from the minimum number of GPS satellites needed 
to resolve the relative position between two antennae. From 
these measurements, a set of all possible integer combina- 
tions is derived. Using carrier phase measurements from an 
additional GPS satellite, the unlikely integer combinations 
are systematically eliminated. This technique is suited to the 
context of attitude determination where both receivers use 
the same frequency standard and the distance between the 
antennae is fixed. This approach, however, is ill-adapted for 
positioning over large displacements, where the initial set of 
satellites is four and the distance between the receivers is not 
known a priori, the technique is then extremely susceptible 
to noise, and computationally intensive. Knight (U.S. Pat. 
No. 5,296,861) details an approach similar to that of Hatch, 
except that a more efficient technique is derived for elimi- 
nating unlikely integer combinations from the feasible set. 
Knight’s technique also assumes that the two receivers are 
on the same clock standard. 

Counselman (U.S. Pat. No. 5,384,574) discloses a tech- 
nique for GPS positioning that does not resolve integer cycle 
ambiguity resolution but rather finds the baseline vector 
between two fixed antennae by searching the space of 
possible baseline vectors. The antennae track the GPS 
satellite signals for a period of roughly 30 minutes. The 
baseline is selected that best accounts for the phase changes 
observed with the motion of the GPS satellites. This 
technique, however, assumes that the baseline vector 
remains constant over the course of all the measurements 
during the 30 minute interval, and is therefore only suitable 
for surveying applications. Moreover, it also assumes that 
the clock offset between user and reference receivers 
remains constant over the 30 minute measurement interval. 

Amotion-based method for aircraft attitude determination 
has been disclosed. This method involves placing antennae 
on the aircraft wings and tail, as well as a reference antenna 
on the fuselage. The integer ambiguities between the anten- 
nae can be rapidly resolved as the changes in aircraft attitude 
alter the antenna geometry relative to the GPS satellite 
locations. This approach, however, is limited to attitude 
determination and is not suitable for precise absolute posi- 
tioning of the aircraft itself. 

Current state-of-the art kinematic carrier phase GPS navi- 
gation systems for absolute positioning have been disclosed. 
These systems achieve rapid resolution of cycle ambiguities 
using ground-based navigational pseudo satellites 
(pseudolites) which transmit either an additional ranging 
signal (Doppler Marker) or a signal in phase with one of the 
satellites (Synchrolites). Although this approach rapidly 
achieves high precision absolute positioning, it provides 
high precision and integrity only when the user moves near 
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the ground-based pseudolites. In addition, the pseudolites 
are expensive to maintain. 

Therefore, each of the existing techniques for satellite- 
based navigation suffers from one or more of the following 
drawbacks: (a) it does not provide centimeter-level accuracy, 5 
(b) it does not quickly resolve integer cycle ambiguities, (c) 
it is not suitable for kinematic applications, (d) it provides 
only attitude information and does not provide absolute 
position information, (e) it does not have high integrity, (f) 
it requires the deployment and maintenance of pseudolities, 10 
(g) its performance is limited to users in a small geographi- 
cal area near pseudolities, or (h) it requires the user receiver 
and/or the reference receiver to have an expensive highly 
stable oscillator. 

15 

III. SUMMARY OF THE INVENTION 

In view of the above, it is an object of the present 
invention to provide a system and method for centimeter- 
level kinematic positioning with rapid acquisition times and 
high integrity. In addition, it is an object of the invention to 
provide such a method that does not depend on additional 
signals transmitted from nearby navigational pseudolite 
transmitter, and does not require a highly stable oscillator 
such as an atomic clock. It is further an object of the ^ 
invention to provide a navigation system requiring only 
carrier phase information. These together with other objects 
and advantages will become apparent in the following 
description. 

In order to obtain high-integrity estimation of integer 30 
cycle ambiguities, carrier-phase measurements must be 
made for a time interval long enough that the displacement 
vectors between the user and the signal sources undergo 
substantial geometric change. Surprisingly, the present 
inventors have discovered a method and system for fast 35 
acquisition, high integrity, kinematic carrier-phase position- 
ing using a non-navigational signals from low earth orbit 
(LEO) satellites which are not necessarily intended for 
navigational use. The short orbital periods of these LEO 
satellites provide the required change in geometry for reso- 40 
lution of cycle ambiguities with high reliability in a few 
minutes. The technique, therefore, provides fast acquisition, 
high precision and high integrity without depending on 
signals from ground-based pseudolites in close proximity to 
the user. In addition, the technique has the advantage that it 45 
does not require that the LEO satellites have any special 
features (e.g. atomic clocks or the ability to transmit navi- 
gational signals). 

Remarkably, the present inventors have discovered a 
system and method for satellite -based navigation using 50 
signals from non-navigational satellites in low earth orbit. 
Beginning only with an estimate of the user clock offset, 
high precision kinematic positioning is provided using only 
carrier-phase signals transmitted from earth orbiting satel- 
lites. By using signals from at least one LEO satellite, high 55 
integrity and fast acquisition is provided. The other signals 
may be from other satellites, including high earth orbit 
navigational satellites, or from any other space-based or 
earth-based sources, including pseudolites and other earth- 
based transmitters. Only the carrier phase signals from these go 
other sources are required. 

In a preferred embodiment, an initial estimate of the user 
position and clock offset is provided by conventional code- 
phase differential GPS techniques. In addition, differential 
carrier phase measurements are used in order to eliminate 65 
errors caused by atmospheric phase distortion, satellite 
ephemeris deviations, and deliberate corruption of the sat- 
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ellite signals. As will become apparent, however, the fun- 
damental technique of using non-navigational carrier signals 
from LEO satellites for resolving integer cycle ambiguities 
in a navigational system is not limited to these specific 
implementations. In alternative embodiments, for example, 
an initial clock offset may be estimated by any combination 
of known techniques for navigation, including anything 
from sophisticated earth-based radio navigation to simply 
calibrating the user receiver to a known reference. 

In the preferred embodiment, centimeter-level positioning 
is provided by combining the navigational data available 
from GPS satellites with the non-navigational carrier phase 
data available from LEO satellites. In addition, the non- 
navigational carrier phase data from the GPS satellites is 
used. The method is robust to frequency-dependent phase- 
lags in the navigation receivers, as well as to instabilities in 
the crystal oscillators of the satellites, and of the receivers. 
We describe how the general technique can be applied to a 
variety of different satellite communication configurations. 
Such configurations include satellite transmission of mul- 
tiple beams with different phase-paths, and bent-pipe archi- 
tectures where the uplink signal is frequency-converted by 
the satellite and retransmitted. 

Generally, in one aspect of the invention a user device is 
provided for satellite -based navigation. The device com- 
prises at least one antenna for coupling to signals transmitted 
from a set of satellites. The set of satellites includes a set of 
LEO satellites that do not necessarily transmit navigational 
information. A receiver in the device tracks the signals to 
obtain carrier phase information comprising geometrically 
diverse carrier phase information from the LEO satellites. A 
microprocessor in the device calculates the precise position 
of the user device based on the carrier phase information and 
an initial estimate of the device clock offset. In a preferred 
embodiment, the device calculates an initial estimate of 
position and clock offset from code phase information 
derived from navigational signals transmitted by naviga- 
tional satellites. In addition, the preferred embodiment uses 
reference carrier phase information transmitted from a ref- 
erence station to improve the accuracy of the position 
estimate. 

In another aspect of the invention, a satellite -based navi- 
gation system is provided. The system comprises a set of 
satellites, including LEO satellites, that transmit carrier 
signals, a reference station, and a user device. The reference 
station samples the carrier signals to obtain reference carrier 
phase information which is then transmitted to the user 
device over a communication link. In addition to receiving 
the reference carrier phase information, the user device 
directly tracks the carrier signals to obtain user carrier phase 
information from the set of LEO satellites. The user device 
then calculates its precise position based on the reference 
carrier phase information and the user carrier information. 
The calculation uses the geometrically diverse reference and 
user carrier information from the LEO satellites to quickly 
resolve parameters related to the integer cycle ambiguities in 
the reference and user carrier phase information. In a pre- 
ferred embodiment, the calculation of the user position is 
based on an initial estimated clock offset and position 
calculated from navigational code phase signals transmitted 
from a set of navigational satellites. Preferably, the reference 
station also transmits differential code phase correction data 
to the user to improve the accuracy of the initial estimate. 

In another aspect of the invention, a method is provided 
for estimating a precise position of a user device in a 
satellite-based navigation system. The method comprises 
transmitting carrier signals from a set of satellites, wherein 



US 6,373,432 B1 


5 

the set of satellites includes a set of LEO satellites; tracking 
at a user device the carrier signals to obtain user carrier 
phase information comprising geometrically diverse user 
carrier phase information from the set of LEO satellites; and 
calculating the precise position of the user device based on 5 
an initial position estimate and the user carrier phase 
information, wherein the calculation uses the geometrically 
diverse user carrier information from the set of LEO satel- 
lites to quickly resolve integer cycle ambiguities in the user 
carrier phase information. In a preferred embodiment, the 10 
method includes tracking at a reference station the carrier 
signals to obtain reference carrier phase information com- 
prising geometrically diverse reference carrier phase infor- 
mation from the set of LEO satellites. The reference carrier 
phase information is then used to improve the accuracy of 15 
the position calculation. In a preferred embodiment, the 
method further comprises estimating an approximate user 
position and clock offset using code phase signals received 
from a set of navigational satellites. Preferably, differential 
code phase techniques are used to improve the accuracy of 20 
the initial estimate. The preferred embodiment of the method 
also includes additional advantageous techniques such as: 
compensating for frequency dependent phase delay differ- 
ences between carrier signals in user and reference receiver 
circuits, reading navigation carrier information and LEO 25 
carrier information within a predetermined time interval 
selected in dependence upon an expected motion of the user 
receiver and the LEO signal sources, calibrating LEO oscil- 
lator instabilities using navigation satellite information, 
compensating for phase disturbances resulting from a bent 30 
pipe LEO communication architecture, compensating for 
oscillator instabilities in the user and reference receivers, 
predicting present reference carrier phase information based 
on past reference carrier phase information, and monitoring 
the integrity of the position calculation. 35 

IV. BRIEF DESCRIPTION OF THE FIGURES 

FIG. 1 shows an operational overview of a preferred 
embodiment of the invention. 

FIG. 2 shows some possible methods of conveying the 40 
LEO ephemerides to the user in a system according to the 
invention. 

FIG. 3 shows a receiver architectural overview according 
to a preferred embodiment of the invention. 

FIGS. 4 a and 4b show two different mixing, filtering and 45 
sampling schemes according to a preferred embodiment of 
the invention. 

FIG. 5 shows tracking and phase counting assemblies for 
the mixing scheme of FIG. 4 a, according to a preferred ^ 
embodiment of the invention. 

FIG. 6 shows tracking and phase counting assemblies for 
the mixing scheme of FIG. 4b, according to a preferred 
embodiment of the invention. 

FIG. 7 shows tracking and phase counting assemblies for 55 
the mixing scheme of FIG. 4b, and a phase latch 
architecture, according to a preferred embodiment of the 
invention. 

FIG. 9 shows a microprocessor block diagram for a user 
receiver according to a preferred embodiment of the inven- 60 
tion. 

FIG. 10 shows a microprocessor block diagram for a 
reference receiver according to a preferred embodiment of 
the invention. 

FIG. 11 is a conceptual illustration of lattice points 65 
defined by a lattice basis G in accordance with a preferred 
embodiment of the invention. 
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FIG. 12 shows the beam arrangement for the Globalstar 
S-band downlink. 

FIG. 13 shows the operational overview of a method 
using Globalstar satellites according to a preferred embodi- 
ment of the invention. 

FIG. 14 is a graph of the fractional availability of the 
Globalstar Constellation above 10° elevation. 

FIG. 15 is a graph of the availability of RAIM for 
protection radii of 110 cm, assuming GPS alone and GPS 
augmented with Globalstar, according to a preferred 
embodiment of the invention. 

FIG. 16 shows a reference receiver and transmitter archi- 
tectural overview in accordance with a preferred embodi- 
ment of the invention. 

FIG. 8 is a block diagram of one channel of a tracking 
module designed for the CDMA signal of equ. (1), according 
to a preferred embodiment of the invention. 

FIG. 17 is a graph of average standard deviations in radial 
position errors for a mobile user, at 5 km and 1 km from the 
reference station, according to a preferred embodiment of 
the invention. 

FIG. 18 is a graph of the evolution of the Probability of 
Selecting the Correct Integer Set, in accordance with a 
preferred embodiment of the invention. 

V. DETAILED DESCRIPTION 

For the purposes of this description, we assume that an 
initial estimate of position and clock offset is derived sat- 
ellite navigational signals and that the navigation satellites 
employed are the Navstar satellites of the GPS constellation. 
It will be appreciated that other present or future navigation 
satellites may be used in other embodiments of the invention 
and other navigation techniques may be used to estimate an 
initial position and clock offset. We assume also the use of 
differential techniques involving a single reference receiver 
and a single user receiver which requires real-time position 
information. It will be apparent, however, that these differ- 
ential techniques are not necessary in order to practice the 
invention. 

FIG. 1 presents a schematic overview of the system. The 
central components of the system are the Navstar GPS 
satellites 1 a-d, the LEO satellites 2a, b, the user receiver 4 
and the reference receiver 3. The user 4 and reference 3 
receivers track the absolute carrier phase of the Navstar 
satellite signals 5a-d together with the absolute carrier phase 
of multiple LEO satellite signals 6a, b. By the term absolute, 
we mean the phase measurement is accumulated over time 
and is not modulus 2 jt. The motion of the LEOS 2a, b causes 
rapid change in the angles — shown in the figure as 0 2 and 
0 2 — between the baseline 7 from user antenna 17 a to 
reference antenna 17d, and the line of site vectors from user 
antenna 17a to satellites 2a, b. The rapid change in the 
line-of-sight vectors enables the user receiver 4 to resolve 
the integer cycle ambiguities on the Navstar satellite signals 
5a-d as well as parameters related to the integer cycle 
ambiguities on the LEO signals 6a, b, and consequently to 
position itself with cm-level precision with respect to the 
reference receiver 3. 

We describe below the steps involved in the general 
technique. It would be clear to one skilled in the art how the 
order of these steps would change, or how some steps would 
be adjusted for a static user, or an attitude determination 
problem where use and reference receivers are driven by a 
common oscillator and are separated by a constant distance. 

The reference 3 and user 4 receiver obtain up-to-date 
satellite ephemeris information. 
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Both reference 3 and user 4 receiver measure the code 
phase delay on the signals 5a-d transmitted by the GPS 
satellites. This measurement is known in the art as raw 
pseudorange. 

Based on the code phase measurement, the user 4 and 
reference 3 correlate their clocks to within 1 /usee of 
GPS time. 

In a preferred embodiment, the reference receiver 3 
calculates differential corrections for the code -phase 
measurements and conveys them to the user receiver 4 
over a communication link 8. The user receiver 4 then 
positions itself with meter-level accuracy relative to the 
reference receiver 3 using the differentially corrected 
code -phase measurements. 

The user 4 and reference 3 receiver simultaneously track 
the absolute carrier phase of GPS satellite signals Sa-d 
and LEO satellite signals 6a, b over an interval of time. 

If necessary, the reference receiver calibrates the LEO 
satellite oscillators using the technique described in 
section (VII.C). 

The reference receiver 3 conveys to the user its carrier 
phase measurements and measurement correction data 
over the communication link 8. 

The user corrects for deterministic errors in the carrier 
phase measurement. 

The user employs the data-reduction technique described 
in section (VLB) to identify the integer cycle ambigu- 
ities on the Navstar satellites signals 5a-d ■, and param- 
eters related to the integer cycle ambiguities on the 
LEOS signals 6a, b. 

Once these parameters are identified, the user receiver 4 
is able to position itself in real time, with centimeter- 
level accuracy relative to the reference receiver 3. 
V.A. Data Links in the System which have Many Possible 
Implementations. 

V.A.l Communicating the LEO Satellite Ephemerides 

The navigation capability of the system hinges on the user 
knowing the location of the LEOS 2a, b to reasonably high 
accuracy — see section (VIII.D). Knowledge of the satellite 
position as a function of time can be obtained via the satellite 
ephemeris data. This ephemeris data consists of several 
parameters for each satellite, describing the satellite orbit 
and changes in that orbit over time. To maintain the desired 
levels of accuracy, the ephemeris data must be updated by 
the user roughly daily. FIG. 2 displays the different mecha- 
nisms by which the ephemeris data may be conveyed to the 
user. Any sequence of arrows leading from the satellite 
operations and control center — SOCC 9, or tracking station 
10 to the user 4 a-c is possible. The position sensing for the 
ephemeris data can be achieved either by position sensors on 
the satellites 2 a-c, such as GPS receivers, or by a tracking 
station 10 processing Doppler information from ground 
receivers at surveyed locations to calculate the LEOS’ 2 a-c 
orbital parameters. Whether the information is attained from 
the SOCC 9, or from a separate tracking station 10, the data 
is conveyed to an ephemeris data provide 11 who must make 
the information accessible to the user 4 a-c. A simple imple- 
mentation connects the reference 3 to the ephemeris pro- 
vider 11 via line modems 12a, c over a regular land telephone 
line 12 b. Alternatively, the information can be obtained by 
the reference over a LEO satellite data link 13 a-b. The 
reference would then convey this information to the user 
4 a-c via a data Link from Reference to User — LRU 8. 
Another embodiment has the LEO satellite 2b broadcasting 
ephemeris data 14 a — e on a dedicated broadcast channel 
14 b-e which is received by both the reference 3 and the 
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users 4 a-c. The techniques for deriving ephemeris data from 
satellite position sensors, or from Doppler-tracking 
satellites, are well understood in the art, as are the imple- 
mentations of these different methods for conveying that 
5 information to the user. 

V.A.2 The data Link from Reference to User (LRU) 

For a mobile user, the LRU 8 is implemented with a 
real-time radio connection 1 Sa-d. For an attitude determi- 
nation problem, the LRU 8 could be implemented with a 
10 real-time cable connection, and for static surveying 
applications, the LRU could be off-line. For the mobile user 
the LRU 8 will be implemented using the same basic 
frequency and modulation scheme of either the GPS signals 
5a-d, or the LEO satellite signals 6a, b. In this way, existing 
15 receiver front-end hardware — see FIG. 3 — can be used for 
the LRU 8. The radio LRU could be implemented by 
transmitting a signal 1 5a, b directly to the user 4a, b via a 
reference station transmitter 16, or by using an existing LEO 
satellite data link 1 5c, d. 

20 The central function of the LRU 8 is to convey carrier 
phase measurements made at the reference station 3 to the 
user 4. If the user receiver 4 knows the location of the 
reference station 3, it is able to accurately predict the phase 
measurements that the reference receiver 3 will make for up 
25 to a few seconds. Consequently, the LRU 8 needs to be 
active every few seconds over the duration of navigation. In 
addition to the carrier phase information, the LRU 8 can 
convey to the user 4: 

The satellite ephemerides. This would be necessary for an 
30 implementation where the ephemeridies are not known 
to the user, and cannot be obtained by a directly 
broadcast satellite data link 14 a-e. 

An estimate of the reference receiver’s clock offset. This 
can be used to correct the differential measurement as 

or 

described in section (VI .A). 

Differential corrections for code-phase measurement. In 
order to achieve an initial differential position estimate 
which is accurate to within meters, the reference 3 
4Q would send to the user 4 a set of corrections for the 
range measurements to improve the code phase perfor- 
mance of GPS. The technique of deriving these cor- 
rections is well understood in the art. 

The position of the reference station antenna 17 d. The 
45 data reduction technique depends upon the user know- 
ing the rough location of the reference station antenna. 
The user receiver 4 finds the centimeter-level position 
between its antenna 11 a-c and the reference receiver 
antenna 17^. Therefore, the universal accuracy of the 
50 user’s derived position depends on the reference 

antenna lid position information. 

Status information on the satellites la-d,2a-c. The ref- 
erence station can also send the user information about 
the health and signal characteristics of each of the 
55 satellites being tracked. 

Error correction information. This information is 
employed by the user to minimize the residual errors of 
differential phase measurement, due for example to 
ionospheric and tropospheric delays, and sever satellite 
60 oscillator instabilities. 

V.B. Description of a combined GPS and LEO receiver 

FIG. 3 illustrates the essential components for a reference 
3 or a user 4 receiver, assuming that transmission are 
receivable from the Navstar Constellation (N) as well as 
65 LEO constellations (L 1? L 2 , L 3 ). Only one LEO constellation 
is necessary for the invention, however 3 LEO constellations 
are assumed for this diagram in order to illustrate the 
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scaleable nature of the invention. While the overall structure 
of this receiver cannot be substantially altered, it should be 
recognized that the individual modules identified in FIG. 3 
can be implemented in a variety of ways. 

The antennae subsystem 17 must be sensitive to the 5 
downlink bands of each constellation being tracked. We 
assume in our description of the data reduction technique 
that the antenna subsystem consists of a single antenna 
which is resonant at the relevant frequencies and has phase 
centers 21 a-d for the respective bands which are separated 10 
by only a few millimeters. This assumption simplifies the 
description, but is not essential to the invention. If an 
antenna subsystem is used with phase centers which are 
substantially separated for the different frequencies, the 
relative position of the phase centers would simply be used 15 
to correct the carrier phase measurements — equ. (13) — to 
take this separation into account. This would involve either 
assuming an orientation of the antenna subsystem 17, or for 
attitude-determination problems, iteratively making the cor- 
rection based on the estimated orientation of the vehicle on 20 
which the antenna subsystem is mounted. 

The antenna feeds a low noise amplifier (LNA) 22 which 
has gain over the full bandwidth of the signals being tracked. 

Of course, separate LNA’s could be used for each of the 
constellations being tracked, namely on the N, LI, L2 and 25 
L3 paths 18, 19a-c respectively. We will examine the 
receiver path for one of the constellations, LI 19a, since all 
paths are similar at this level of description. The signal from 
the LNA 22 is bandpass filtered with an rf bandpass filter 23 
and then downconverted by mixing with a locally generated 30 
rf frequency f r $24. An if bandpass filter 256 at the mixer 
26 b output removes the unwanted upper sideband. A 
microprocessor-controlled automatic gain control 29 then 
adjusts the magnitude of the signal to achieve optional use 
of the available sampling bits. For a situation where a signal 35 
is hidden in bandpass noise substantially larger power — 
such as arises with spread spectrum modulation methods — 
the SNR forfeited by sampling with one bit is roughly 1.96 
dB. Therefore, it is possible to implement this system with 
only one bit, in which case the controllable gain of the AGC 40 
is unnecessary. Another potential variation on this architec- 
ture is to place a low noise AGC at the input to the rf 
bandpass filter 23,30-32 for each of the receiver paths, to 
substitute for the LNA 22 and AGC’s 29,33-35. 

After amplification by the AGC 29, the signal 53 enters an 45 
if mixing stage 27, followed by a filtering stage 36 and a 
sampling stage 37. The signal is mixed in an if mixer 27 with 
a locally generated /^frequency f £f 28. The exact means of 
downconverting, filtering and sampling the signal at the if 
stage varies substantially with the signal structure and 50 
designer’s preference. FIGS. 4 a and 4b illustrate two pos- 
sible means of implementing the if mixing, filtering and 
sampling. Both of these schemes would be possible for 
Binary-Phase-Shift-Keyed (BPSK) and Quadrature-Phase- 
Shift-Keyed (QPSK) signals. In the scheme of FIG. 4a, only 55 
a single mixer 40 and filter 41 are used, and the sampling 
section 42 outputs a single sample 43. If the incoming signal 
53 had quadrature modulation, FIG. 5 describes a possible 
architecture for a tracking assembly 44 b preceded by the 
scheme of FIG. 4a. The tracking modules 48a-c could mix 60 
the incoming signal 43 with both an in-phase and quadrature 
component of the output of the Numerically Controlled 
Oscillators 49 a-c, in order to isolate the in-phase and 
quadrature modulation. In the scheme of FIG. 4b, both an 
in-phase 45a and quadrature 45b mixer are used, each of 65 
which outputs to a separate filter 4 6a, b and sample 47 a,b. If 
the incoming signal 53 has quadrature modulation, FIG. 6 


describes a possible architecture for a tracking assembly 44 b 
preceded by the scheme of FIG. 4b. The tracking modules 
52a-c of the tracking assembly could mix the incoming I 
50a and Q 50 b signals with a single in-phase output of the 
Numerically Controlled Oscillator 51a-c to isolate the 
in-phase and quadrature modulation. 

The bandwidth of the filters 41, 46a, b is chosen to 
accommodate the full bandwidth of the signal tracked, 
which is shifted by the offset f 0 =f L 1 S r$ L1 ~^ if L1 ■ These filters 
41, 4 6a, b could be lowpass or bandpass, depending on t 0 . 
The sampling rate should be roughly 5 to 10 times the 
highest frequency component in the signal and the number 
of bits sampled could vary from 1 to 16 bits, depending on 
the signal structure, the SNR, the desired robustness to 
interference, and hardware costs. 

Since the tracking assemblies of FIG. 5 and 6 have 
essentially similar structure, we will consider in more detail 
only the architecture displayed in FIG. 6. We assume the if 
mixing is performed using the scheme of FIG. 4b. The thick 
black lines represent digitized 1 50a and Q 506 samples. The 
incoming samples are latched 12a, 12b and input to S 
tracking modules, where S is the maximum number of 
signals from a particular satellite constellation that one seeks 
to track. Each tracking module 52a^c tracks one satellite 
downlink signal by means of a phase-locked loop. Many 
different techniques for implementing the tracking modules 
52a-c are known in the art. Whichever tracking module 
architecture is employed, it involves an oscillator which is 
phase-locked to the phase of the incident signal. The pre- 
ferred embodiment employs a numerically controlled oscil- 
lator 51a^c in the phase locked loop. In order to analyze 
phase-tracking behavior, we present one possible design of 
a tracking module for generic CDMA signal of the form 


s(r) = — — D(r)C/(r)cos(ti>r + <p{t)) 4 
V2 


—j= D{t)C Q (r)sin(wr + <p{tj) + n(t ) 


( 1 ) 


where D(t) refers to the outer data sequence modulated on 
both the in-phase and quadrature signals. C/t) and C e (t) are 
respectively the in-phase and quadrature spreading 
sequences. n(t) Represents thermal input noise, which is 
assumed to be normally distributed, of zero mean and of 
uniform spectral density N 0 . A tracking module designed for 
the signal structure of equ. (1) is illustrated in FIG. 8. 
Ignoring the effort of front-end gains equally applied to 
signal and noise, the in-phase 98a and quadrature 986 digital 
signal entering the tracking module can be described: 

A (2) 

h = —^=D k Ci k cos(<p k ) + l nk 
V2 

A 

Qk = -^==D k C Qk cos(<f> k ) + Q nk 
£<&) = £{<&>= ^ 

where B c is the pre -correlation signal bandwidth, deter- 
mined by the filters 46a, 6. The upper sideband emerging 
from the carrier mixers 100a, 6 has frequency ~2f 0 and is 
rejected by the accumulator lOla-d which has effective 
bandwith 1/T, where T is the period of the inner codes C/t) 
and C G (t). Consequently, we consider only the lower side- 
band of the mixer outputs 85a, 6: 
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hk — C/kCOS^fik — <p r k) + hnk 

V2 

A 


Qlk — 7=Dk CQkCOS^k — (fork) 4- Qlnk 
V2 


(3) <4t = E{ll(g i }-E l {l2 i Q2i} 

_ K*No, No _ x 
r A 2 T V + A 2 t) 


(7) 


EVL} = E{Qik}' 


No B c 


Each of these signals is then mixed with a prompt inner 
code 102 a,b and a tracking inner code 103 a,b, which con- 
sists of the difference between early and late code replicas, 
separated by some number of chips, d, where d<2. The 
prompt accumulator 101fr,c outputs can be described: 


By monitoring the signal amplitudes of the correlator 
outputs 104 a,b, and varying the loop control accordingly, a 
specific phase-lock loop transfer function, H, is maintained 
10 by the microprocessor 56. Generally, in selecting the band- 
width of the loop, a tradeoff must be made between rejecting 
thermal noise on the one hand, and tracking performance on 
the other. For loop transfer function H, the NCO phase error, 
which we cannot directly estimate by the technique 
15 described in section (VI.B), has variance 


A 

h = COS (<p k - (prk) + l2ni 

V2 jt =1 


°% ~ \ H Uu)\ 2 du+ s o ((o)\l- H(ja ))\ 2 du (8) 


Qli = -j=R{Ti)D^ Sin(^ -<p n t) + Q2ni 
V2 k=l 




MN 0 

~xf7 


where S n (0) is the power spectral density near the origin of 
the thermal noise. The expression assumes that the loop 
bandwidth, 


where R(t,) is the cross correlation between the incoming 
code and the generated code for a time misalignment of x t . 
R(t)» 1 once code-lock is achieved . For an inner code 
period, T, and sample rate f s , the summation is typically over 
M=T% samples. This number may be varied to accommo- 
date code Doppler. Assume that over one correlation period, 

T, the NCO 105 suffers a constant frequency error, Af I -«l/T, 
and a constant phase error, A<j> f -. We can then treat the 
summation in equ. (4) as a continuous integral, and find 

E{I 2i } « — — DiSmciTjiAfiTfcosiAdi) 

V2 

E{Q 2i \ ~ —= D t sinc(27r A/j- r)sin( A0/ ) 

V2 

For an assumed orbital altitude of 1400 km and transmis- 
sions in the S-band, we expect a maximum phase accelera- 
tion due to Doppler of ~100.2 jt rad/s 2 . For T=1 ms, we 45 
expect sinc(2jtAQ>0.97 so the factor can be safely dropped. 
The samples 104 a,b are input to the microprocessor 56 
which estimates the phase error and implements a loop filter 
to achieve desired phase-lock loop performance. Since the 
signal of equ. (1) has a common outer data sequence 86 on 50 
the in-phase and quadrature components, a simple Costas 
Loop discriminator approximates the phase error by multi- 
plying the samples 104 a,b. 

8<pi = hiQn < 6 ) 55 

a 2 m 2 

Em) = A <Pi 

From equ. (6), we see that the gain of the discriminator is 60 
a 2 m 2 


In addition, the variance of the discriminator output can be 
computed. 


(5) 

40 


= I HU<o)\ 2 da>, 

is much smaller than B c . Since the correlator bandwidth, 

1 

T' 

is also much less than B c , we can safely approximate 
S n (0)=Ta2/6c|). From equ. (7), the phase error variance in 
equ. (8) becomes 

y 2 N 0 B l / N 0 \ (9) 


This thermal noise variance is not heavily dependent on 
there being the same data sequence, {D-J, on in-phase and 
quadrature signal components. For example, for a CDMA 
signal with QPSK data modulation, we would use a fourth- 
power loop, rather than the Costas Loop. It has been shown 
that thermal noise variance could then be approximated 


2 N 0 B l 
A 2 



A phase-counter 57 a-c keeps track of the absolute phase 
of the NCO 51 a-c locked to each signal. The phase mea- 
surement of the phase counter SI a-c contains both a frac- 
tional and an integer component, where each integer refers 
to a full 2 jt phase cycle. There are multiple different means 
by which the microprocessor 56 can read the phase of the 
phase counters 57 a-c. Whichever method is employed, for 
each reading epoch, the time from reading the first phase 
counter tracking a signal to the last phase counter tracking 
a signal should span an interval less than a few //sec in order 
to cause errors below the experimental noise floor. This 
specification is made concrete in section (VI. A). One 
method of satisfying this specification is illustrated in FIG. 
6. The microprocessor 56 sequentially reads the phase of 
each tracking module that is locked on a signal. This is 
achieved by means of a select signal 58 which is input to the 
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bus interface 59 of each of the phase counting assemblies 
53 a-d. Each select signal selects the output of one of the 
phase counters in one of the phase counting assemblies. The 
clock driving the phase 86 and select 58 buses should be fast 
enough to enable reading of all of the relevant phase 
counters in all of the phase counting assemblies in the 
specified interval. Another approach to the problem is illus- 
trated in FIG. 7. In this configuration, a single latch signal 
is used to latch 64a-c the phase of each of the counters 
62 a-c in all of the phase-counting assemblies 53 a-d simul- 
taneously. The latched data 61 a-c can then be read sequen- 
tially over a longer interval, by sending select signals 58 to 
the bus interface 63 for each of the relevant latches. Phase 
reading epochs, in which the phase data for all of the signals 
being tracked is read, occur at frequencies roughly between 
1 and 100 Hz, depending on the requirements of the specific 
application. 

The frequency synthesizer subassembly 65 for the 
receiver 3,4 is driven by a single crystal oscillator 66. Many 
different techniques are known in the art for generating a 
series of frequencies from the crystal oscillator input 67. 
One method of generating the desired frequencies of a 
receiver path for the Globalstar LEO constellation is dis- 
cussed in section (VIII). Areal-time clock 68 is maintained 
continuously, and generates time information for the micro- 
processor 56 to aid in the initial acquisition of satellites. 

The microprocessor functions and memory for the user 4 
and reference 3 receivers can be decomposed as illustrated 
in the block diagrams of FIGS. 9 and 10. Each micropro- 
cessor consists of a CPU 73a, 5, memory 70a, 571a, 5 and a 
bus interface 74a, b through which to communicate with the 
rest of the receiver. The memory of the microprocessor 
stores the software routines 70a, b as well as the data 71a, b 
necessary for the implementation of those routines. We will 
first consider the software routines 70a for the user receiver. 

70a. 1 — The initial acquisition algorithm performs a sig- 
nal search routine to initially establish phase-lock on satel- 
lite signals. This involves using the satellite ephemeris data 
71a. 1 , the GPS signal structure data 71a.2, the LEO signal 
structure data 71a. 3, possibly the location of ground uplink 
stations 71a. 4, together with data 75a-d from the tracking 
assemblies 44a-d to control the tracking modules such that 
phase lock is established on the available satellite signals. 
Control commands are applied via the select/control bus 76. 
The initial acquisition can also involve implementing 
frequency -locked loops as well as delay -locked loops, and is 
well -understood in the art. 

70a.2 — The maintenance of phase-locked loops involves 
controlling the components of the tracking assemblies 44a-d 
so that phase lock is maintained on all relevant signals. The 
microprocessor implements a control law for closing the 
phase-locked loops, and possibly, delay-locked loops for 
each of the tracking modules in use. These control tech- 
niques are well-understood in the art. 

70a. 3 — Interpreting and demodulating data from the sat- 
ellite downlinks involves using the data about the LEO 
71a.3 and GPS signal structures 71a.2 to read and interpret 
the outputs of the tracking modules 7 5a-d. 

70a. 4 — The carrier phase measurement routine is the 
process whereby carrier phase data is read from the phase 
counters of each of the phase -counting assemblies 53 a-d 
and interpreted to produce phase measurements that can be 
input into the data-reduction routines 70a. 6. 

70a. 5 — The code phase measuring and positioning rou- 
tine is the process whereby code phase data is read from the 
GPS tracking assembly 44a, corrected using the differential 
corrections received from the reference and processed to 
obtain a meter-level position estimate and clock offset 
estimate. 
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70a. 6 — The position calculation routine involves the cor- 
rection of carrier phase measurements as well as an estima- 
tion algorithm. The correction of carrier-phase measure- 
ments is the process whereby carrier-phase data is corrected 
5 for deterministic disturbances based on information received 
from the reference, satellite ephemeris information 71a. 1, 
possibly the location of ground uplink stations 71a. 4 and 
possibly error prediction data 71a. 5. The nature of these 
corrections are discussed in sections (VI.A, VILA, VII.B). 
10 The estimation algorithm is the key process by which data 
is processed to identify the integer cycle ambiguities for the 
GPS satellites as well as related parameters for the LEO 
satellites, and to subsequently position the user with 
centimeter-level precision relative to the reference receiver 
15 3. The algorithm employs the satellite ephemeris informa- 
tion 71a. 1, the GPS signal structure data 71a.2 the LEO 
signal structure data 71a. 3, and possibly the location of 
ground uplink stations 71a. 4. The algorithm is discussed in 
detail in the section VLB. 

20 70a. 7 — Receiver autonomous integrity monitoring may 

be used by the receiver to independently check the validity 
of the position solution using the satellite ephemeris infor- 
mation. This technique is well-understood in the art and is 
described in more detail in section (VIII.B). 

25 70a. 8 — A phase velocity measurement routine may be 

employed to enhance the position estimation methodology. 

70a. 9 — Filtering position estimates involves applying to 
the position data a digital filter, which takes into account 
known aspects of the users’ motion, such as bandwidth 
30 constraints, in order to generate more accurate, noise-free, 
position estimates. This could also involve the use of 
Kalman-filtering techniques to combine the carrier-phase 
position estimates with data obtained from other sensors 
such as accelerometers and gyros. 

or 

Most of the routines discussed also makes use of assorted 
miscellaneous variables 71a. 6 We turn out attention now to 
the reference receiver, to discuss those routines which are 
not necessarily applicable for the user receiver. 

70/). 5 — The LEO clock calibration routine is used to 

40 

identify frequency offsets of the LEO satellite oscillators. 
The algorithm makes use of the satellite ephemeris infor- 
mation 715.1, the GPS signal structure data 715.2, the LEO 
signal structure data 715.3 and possibly the location of 
45 ground uplink transmitters 715.4 to identify the frequency 
offset of the LEO downlink due to the long-term instability 
of the satellite oscillator. The technique is detailed in section 
(VII.C). 

705.6 — The code phase measurements and differential 
50 correction calculation involve reading the code phase mea- 
surements from the BPS tracking assembly 44a and com- 
bining this information with the known position data 715.6 
for the reference antenna 17 d, to calculate differential cor- 
rections for the user. 

55 705.7 — Transmitting data to the user is the process of 

coding and transmitting the data destined for the user 
receiver 4 in terms of the data communications protocol of 
the LRU 8 that is employed. 

We have described the microprocessor operation assum- 
60 ing that all of the code is implemented on a single micro- 
processor. Another implementation might have multiple 
microprocessors in the receiver, each with specialized tasks. 
For example, a microprocessor for each tracking assembly 
might maintain signal tracking loops, while a separate 
65 microprocessor would be dedicated to the computation- 
intensive data reduction routines. Multiple permutation on 
this theme are possible. 
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VI. PRECISE POSITION CALCULATION 
VI .A. Details of the Differential Carrier Phase Measurement 
The scenarios in which centimeter-level positioning accu- 
racy is required within a few minutes of system initialization 
includes surveying, construction, precise control of land 5 
vehicles, as well as high-integrity tasks such as attitude 
determination and automatic landing of aircraft. The data 
reduction technique is similar in all cases, with simplifica- 
tions for the context of a user receiver that is stationary in 
Earth-Based-Earth-Fixed (EBEF) coordinates. We will 10 
focus our attention on the more general case of a mobile user 
receiver 4, which is driven by a oscillator 66 that is distinct 
from that of the reference receiver 3. The purpose of this 
section is to characterize the software and hardware upon 
which centimeter-level navigation accuracy depends. 15 

We denote the nominal carrier frequency of the satellite 
downlink as co 5 . The phase of the satellite frequency syn- 
thesizer at time t s can be described: 

i|i s (f s )=(o s t s +T s A(o s ( a )^ a (10) 20 

where Aco s (Q models the deviation from nominal frequency 
due to drift of the crystal oscillator onboard the satellite. We 
model this drift as a time-varying offset in the satellite clock 

1 r*s 

T s (t s ) = — I A oj{a)da. 


Consequently, 

Similar clock offsets occur at the user and reference 
receivers. At true times t u and t T the user and reference 
receivers respectively record times t u +x u (t u ) and L+x^L). 
Therefore, the phase output from the crystal oscillator 66 of 35 
the user’s receiver 4 at time t u is 

( 12 ) 

where co^ is the nominal oscillator frequency. The frequency 4Q 
synthesizer 65 of the receiver generates the rf and if mixing 
signals by multiplying ikx(t M ) by factors a rf and a if respec- 
tively. We denote the phase of the satellite signal emerging 
from the user receiver’s LNA 22 at time t u as l P s:M (t M ). The 
output of the first mixer 2 6a-d, after bandpass filtering 
25 a-d, has phase 'I , 1 (t„)=T , sl ,(t„)-a Tf m^[t„+t„(t„)]. We 
assume for simplicity of explanation that the scheme of FIG. 

4 a is employed. Hence, the second mixer 40, after filtering 
41, generates a signal with phase T f 2 (t M )= 1 P SM (t M )-(a T:r +a I - r ) 
co^[t M +T M (t M )]. This phase is tracked by the phase-locked loop 5Q 
of the tracking module for satellite s, and is read from the 
corresponding phase counter by the microprocessor 56. 
Since the nominal satellite frequency is co s , the nominal 
offset frequency of the signal tracked by the PLL is co 0 =co s - 
(a^+a^ooy. The microprocessor 56 differences the phase it 55 
reads, «)> w ith the phase component caused by this 
offset frequency, co 0 . Since the microprocessor’s measure of 
time interval is directly affected by crystal 66 instability, it 
calculates this phase component as w 0 [t M +x M (t M )]. The result- 
ing phase measurement is 

<Psu{t u ) = ^0 + T«(r«)] - ¥ 2 (r„) - N su 2tt 

= (x) s \t u + — 'Psh(?h) — N su 2t: 

where we have included an integer cycle phase ambiguity, 
N su , since the microprocessor’s initial phase measurement is 


60 

( 13 ) 

65 


25 


( 11 ) 


30 


16 

modulus 2 jt. Consider now the phase of the incident satellite 
signal, The signal is affected by phase disturbances on 
the satellite-to-user path, as well as frequency-dependent 
phase lags in the receiver. In addition, the phase measured 
depends on the position of the satellite, r s , at the time of 
transmission, rather than at the time of reception. Applying 
these factors to equ. (11), the phase of the signal from 
satellite s emerging from the LNA 22 of the user receiver 4 
at time t u is: 


/ \ I 1 ( Psu{t u ■ 

%u(t u ) = wJ t u - -p su [t u 


( 14 ) 




■ Psu - n su 


where /u su is the frequency-dependent phase delay of the 
receiver and v] su is the error due to ionospheric and tropo- 
spheric delay as well as thermal noise and imperfect carrier- 
phase tracking in the receiver. r\ su Is actually time-varying, 
but we ignore this time dependence for now. In order to 
represent the signal path length, we have denoted by (t 0 ) 
the distance from user’s current position at time t u to the 
satellite position at transmission time t 0 , p SM (t 0 )=|r M (t M )-r s 
(t 0 )|. Since we cannot know exactly the location from where 
a satellite transmitted if we do not know t 0 , the precise 
calculation of t 0 requires an infinite regression. However, 
one can make the simplifying approximation. 



C 


1 

tu- -Psu(tu) 


( 15 ) 


which generates a worst-case ranging error <2 mm for 
satellites at 1400 km (the Global-star nominal orbital 
altitude). Therefore, from equ.’s (13) and (14) we estimate 
the user’s phase measurement at time t M . 


<Psu{tu) = 


- PsuUu - 


Psu(t u )\ 


. . / Psu(tu) y 

w s |T„(r„) - T s \t u J 


( 16 ) 


N su ln + p su + n su 


The receiver assigns a timetag t to the measurement made 
at time t u . Since the user receiver does not know true time 
t M , the timetag will be effected by the clock offset of the 
receiver. The measurement must therefore be recast in terms 
of the clock of the user’s receiver. If a user receiver makes 
code-phase measurements on the GPS signals, the clock 
offset in the receiver, x u (t u ), can be estimated to within 1 
/usee. Once x M (t M ) has been estimated, two algorithmic 
approaches are possible. Firstly, one can use this estimate to 
select the times at which phase data is read from the phase 
counters in the receiver, so as to continually correct for the 
receiver’s clock offset. This clock steering technique limits 
the magnitude of the timetag error, |t-tj, to roughly 1 /usqc. 
Secondly, one can use the estimate of x u (t u ) to actively 
correct for any errors which would arise in the differential 
phase measurement. Since these approaches are conceptu- 
ally very similar, we will describe the latter approach in 
detail. Once the issues are identified and resolved, it will be 
clear to one skilled in the art how the algorithm varies for the 
former approach. 

Assume the user receiver 4 estimates its clock offset to be 
x u using code-phase measurements. We assume for now that 
x u is time -independent since it need not be continually 
updated. We define Ax u (t u )^z u (t u )-x u . For phase data read at 
true time t M , the user receiver assigns a timetag t=t w +Ax w 
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(t M )+T M , where T u is a residual error in sampling time, 
distinct from the clock bias, resulting from the digital logic’s 
imperfect precision in reading phase at the particular instant 
of time identified by the microprocessor. Recasting the 
measurement in terms of timetag t: 


We may rewrite the terms involving 


dpsu 

~dF 


(x) s ( 1 A 

<psu(tu) = — Psul t - A r„(r„) - T u - -p su (t - A r u (t u ) - T u ) I + 
u s r u (t - A r„(r„) - T u ) - N su 2n + p su + n su - 

- Ar„(r„) -T u - ^ p su (t - A r u (t u ) - T„)j 


(17) 


and 


10 


dpsu 

~dF 


as 


Since |Ax M (t M )+Tj is of the order of a few /usee, we may 
Taylor-expand to first order in Ax M (t M )+T M , with negligible 15 
error in dropping the higher-order terms, 


<Psu(tu) ' 


U s ( Psu{t)\ 
— Psu[t + 


( 0 - 


u s dpsuii) 
c dt 


[A r u (t) + T u \- 


U S T S U 


Psu(t) \ 


^Vs« + Psu 3" tl S u 


(18) 


20 


U s dPsr(t ) 

c dt 


\j r (t) - T u (t)] - 


u s dp sr (t) 
c dt 


[tr-? u ] + 


(21) 


u s dp sr (t ) fdp sr (t) 

C dt Vr mJ ( dt 


dPsui* ) 
dt 


j-^-[Ar H (r) + T u \ 


or equivalently as 


u s dp su (t ) 

Mb 

c dt 


t«(0 ] - 


U s dp su (t ) 
c dt 


[t, - T„] + 


( 22 ) 


U s dp su (t) f d Psr (t) dp su m u s 

Note that for highly unstable oscillators, one might 25 c dt 1 r ui { dt dt ) c A r 

include in this expansion the terms 


dr u (t ) 
5 dt 


[A r u (t) + T u ] 


30 


and 


dr s (t ) 
5 dt 


[A T K (t) + T H ], 


35 


which could be incorporated into the estimation algorithm. 
However, these terms are negligible in most implementa- 
tions of the invention. 

A completely similar approach to the reference receiver’s 40 
phase measurement yields the expansion 


, , . U s / Psr(t ) \ . . 

<Psr(t r ) ~ — Psr[t + OJ s T r (t) ' 

c v c > 

U s d p sr (t) 


(19) 


c dt 

Psr + n sr 


-[Ar r (r) + r r ]-w s rJf 


Psr(t) \ 


N sr 2jt -t 


45 


The reference receiver’s phase measurement is time- 50 
tagged with a time t and transmitted (or otherwise 
communicated) to the user. The user matches the timetags on 
the data and performs a single difference, which we now 
index with the timetag t rather than a true time. The 
differential measurement is then 55 


<Ps(t) = <Psuih) - <p sr (t r ) 


( 20 ) 


$ / Psu(t)\ U s / Psr(t)\ 

- Psu {t-—)--Psr{t-—)l 


- r s (r - ^ )] + w s [r„(r) - r r (0] - 

u s dp su (t ) u s dp sr (t) 

— [Ar„(r) + T u ] h — — [Ar r (r) + T r \ - 

c dt c dt 

(N su ~ N sr )2n + (p su ~ Psr) + (n su ~ n sr ) 


60 
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We will adopt the representation shown in equ. (21). It 
will be clear to one skilled in the art how the issues we 
describe can be transferred to equ. (22). The 1 st term in equ. 
(21) can lead to large errors and is incorporated into the 
estimation strategy. The 2 nd term in equ. (21) can be directly 
calculated, and subtracted from the differential measure- 
ment. If the user and reference receivers are implemented 
with similar digital logic, the term |T t -TJ can be made less 
than 0.1 /usee . Hence the 3 rd term can be ignored with 
distance -equivalent errors <1 mm. Now consider the 4 th 
term. For a satellite at 1400 km, and a stationary user 
receiver 10 km from the reference, the term 


dPsrit) 

dt 


dp su (t)\ 
dt ) 


< 50 m/s. 


<50 m/s. To ensure that the 4 th term produces a worst-case 
distance -equivalent error <1 mm we must have |Ax M (tJ+ 
Tj<20 juscc, so that the term can be safely ignored. Consider 
the phase-reading scheme displayed in FIG. 6. The time 
interval required to read all of the active phase -counters of 
all of the active phase counter assemblies 53 a-d contributes 
to the magnitude of the differentially uncompensated term 
|Ax M (t M )+Tj. Hence, for a static user, this time interval 
should be roughly <18 //sec. However, for a mobile user 
receiver moving at, say, 250 m/s, the term 


dpsr(t) 

dt 


dt ) 


< 300 m/s, 


<300 m/s, and the time period should be roughly <3 //sec. 
That rate at which the estimate x u needs to be updated so that 
Ax M (t) remains small depends on the stability of the receiver 
oscillator 66. For example, for a long-term oscillator stabil- 
ity of 1:10 7 , updates every 2 minutes and every 20 seconds 
are sufficient for the static user and mobile user respectively. 

The terms in equ. (20) relating to the satellite clock offset 
can be expanded to first order 
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r ( Psr(t)\ / Psu(t)y\ dr s (t)0) 

— JJ * -^-(PM-Psuit)) 


and we redefine the parameters we seek to estimate as 
(23) follows: 


For a satellite oscillator with long-term frequency stability 
of 1:10 6 we expect 




to get as large as 10 -6 . This could cause distance -equivalent 
errors as large as 1 cm. In section (VII. C), we describe a 
technique for calibrating the frequency offset of the LEO 
satellite oscillator so that the first-order expansions of equ. 15 
(23). can be calculated and subtracted out of the phase 
measurement. 

Eliminating all terms from measurement equ. (20) which 
are either negligible, or actively subtracted out of the 
measurement, the estimate of the resulting measurement is 20 


^s / p su (t)' (x) s / p sr (t)\ 
's(t) ~ Psu[t pjt 


4 1 " ? - Tr(t)] - (Nsu - Nsr)2jr + 

(Psu - Psr) + (n S u ~ n sr ) 


For the Navstar satellites, the new parameter N s simply 
reduces to N s -N a since all GPS satellites are transmitting 
similar frequencies. For the LEO satellites, the phase-delays 
do not cancel the parameters N s cannot be regarded as 
having integer values. If 


| Ni - Ni - — | < 200 cycles, 


<200 cycles, the redefinition (31) leads to a maximum 
distance -equivalent error of <1 mm for satellites at 1400 km. 
The measurement can then be approximated: 


/ Psu(t)\ ( Psr(t)\ 
y s (t) = PsuU — - Psr[t — 


{ c dt ) 


\cT(t)-N s X s +n s 


To clarify the estimation strategy, we redefine the com- 
ponents of this measurement as follows: 


where N ± is by definition 0. 

The set of time-dependent parameters which we seek to 
estimate is 


N=N-N„ 


( 25 ) ®{t)=Vr u {t) T cx{t)Y (33) 

(26) We create an observation matrix for these parameters 

(27) 35 b ase d on our estimates of the line-of -sight vectors to the 

satellites 


n s = j^(n su -n sr ) 


where X s is the nominal wavelength of the satellite carrier. 
Multiplying equation (24) by 


\ L Psu(t)) T 

M c \ t 1 

. PJO ) c 


c dt 


For a current set of estimates, 0(t) and ~N S , we can 
construct an estimate of our prediction error for satellite s: 


to convert phase to a distance measurement, we then have 


.. ( Psu(t)\ ( Psr(t)\ 

y s (t) = p S u[t PsrU 


(, 1 dPsu(t)) . , .. . 

[ 1 -c-arr (,) ~ NsXs + ^ s+ns 


VLB. Estimation Strategy 

If one attempts to estimate all of the integer components, 
{N s }, as well as user position, r M (t) and clock biases x(t), the 
complete set of parameters would be almost unobservable. 
The resulting estimation matrix would be poorly conditioned 
and highly susceptible to measurement noise n s . 
Consequently, we select one of the Navstar satellites say 
satellite 1, to be a reference satellite for differencing. We 
make an initial approximation of the associated integer, N 1 
using equ. (24), based on our estimate of position and clock 
offset using code -phase measurements. Then, we adjust the 
measurement. 


■4-¥)- 


CT{t) - ~N S X S - y s {t) 


Estimation matrices and prediction errors for all visible 
satellites are stacked into combined matrices, 



h{t) 


A ^i(r) ' 


hit) 


A y 2 (t) 

m = 

. hit) _ 

m> = 

.A y s it)_ 


where S is the total number of satellites visible. The batch 
measurement equation, involving measurements from t 1 
through t^ and the batch parameter update, A0, can then be 
expressed as 


ys(i)=y s (t)+fii 



21 


US 6,373,432 B1 


22 


with the matrix structures 


■ Ay(/i) 


A<&( fl )' 

A©(f 2 ) 

A Y(t 2 ) 

A0 = 

A 0(f/v) 

.A Y{t N )_ 


A N 2 



. A N s 


H{h) 0 ■■■ -A 

0 H(t 2 ) 0 ■■■ -A 

0 ■■■ H(t N ) -A 


(38) 


5 


10 


Hdt): 


4-¥) T 

4-“) 



HlU) = 


1 dPlaW 
c dt 


1 dp SuW 

c dt 


( 41 ) 


15 


20 


We then structure the batch 

Hiih) H 2 (h ) 
H Y {t 2 )G{t 2 - h ) 0 

H = 

Hi(t N )G(t N -tO 0 


estimation matrix, H, as 

0 -A' 

H 2 (h) 0 ••• -A 


(42) 


0 


H 2 (t N ) -A 


where 


and the batch parameter update matrix, A0 as 


- 0 


... 0 ' 

a 2 

0 

... 0 

0 

A 3 

... 0 

0 


0 

. 0 


... 


The disturbance matrix V contains errors due to each 
satellite’s measurement n s (t) — which we may reasonably 35 
assume is uncorrelated with distribution n s (t)~N(0,a^ 2 ) — as 
well as ephemeris errors e s (t) due to imperfect knowledge of 
the satellite’s position which affects calculation of p SM and 
p^. Combining these two noise sources, the disturbance 
matrix has the form: 4 


A©=[Ar c (f 1 )Ar y (f 1 )Ar z (f 1 )cAT(f 1 ) . . . cAr(t N )A~N 2 . . . A ~N S Y (43) 

and we may proceed with a batch measurement equation as 
in (37) above. 

Once the estimation problem has been formulated in the 
manner of equ. (37) the solution can be well-conditioned. 
The good conditioning is due to the geometric diversity 
resulting from the motion of the LEO satellites. This geo- 
metric diversity decreases the condition number of the 
estimation matrix H 


k(H) = 


O’max (H) 
O’min (70 


(44) 


«i hi) + ei(ri) 


n s (h) + e s (t 2 ) 


V = 


«i(r 2 ) + ei(r 2 ) 


n s (t 2 ) + e s (t 2 ) 


where a miM (H) and o max ( H) are the minimum and maximum 
(40) singular values of H. The condition number indicates the 
sensitivity of the parameter solution to disturbances V as 
45 well as to errors in the estimation matrix, 8H. This concept 
can be made more concrete by considering the || || 2 norm of 
the parameter estimate errors for a simple least-squares 
parameter solution. Imagine A0* is the parameter update 
solution of the least-squares problem with these error 
50 sources removed: 


It should be noted that the matrix structures and param- 
eters described above can be altered if the user receiver 4 is 
static relative to the reference receiver 3. This situation 
pertains, for example, in surveying applications, or any 
problem where a vehicle can remain stationary until good 
integer estimates are available. In such scenarios we need 
only estimate the 3 coordinates of r M (t 1 ). Given an estimate, 
fj/T) of the initial position, our estimate of the position at 
time t n is simply G(t M -t 1 )r M (t 1 ), where G is a rotation about 
the z-azis in EBEF coordinates, which accounts for the 
earth’s rotation between time t n and t 1 . To account for the 
reduction in the number of parameters, we define for each 
time t two separate stacked observation matrices: 


A0* = arg min \\(H - (5//)A0 - (A Y - V)\\ 2 

AQg ^ 4N+ s ~1 


55 

while A0 i5 is the actual least-squares solution found. 
A0 L * = arg min ||//A0 - Ay|| 2 

A©eR 4/V+5,_1 

60 


It can be shown (see Golub [18]) that if 


65 


\m\ 2 imi 2 ) 

11//M2 ’ iiahi 2 i 


<k(H) 


(45) 


(46) 


( 47 ) 



and 
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. || H&Ols - AK|| 2 

S11 ^ = IIAMI, 


then 


(48) 


5 


E{e s (h)e s {t 2 )} 


^2 t' su {t 2 ) 

e [ r su {h)r su {t 2 ) 


r sr (h) T r sr {t 2 ) 
r S r{h)r sr {t 2 ) 


fsuihf f sr (h) _ r sr (h) T r su (t 2 ) ' 
rsuih )r sr (t 2 ) r sr {h )r su (r 2 ) 


= o 2 es {t i , t 2 ) 


( 55 ) 


||A0 LS -A0J| 2 

HA0ls||2 




(49) 

10 


For one skilled in the art, equ. (49) indicates how small 
k(H) should be for a given e in order to achieve a desired 
level of accuracy in the parameter estimates. 

VI.C Mathematical methods for solving the estimation prob- 
lem 


Many different mathematical methods can yield a solution 
to the problem posed in equ. (37). A method can be chosen 
depending on the processing power available in the receiver 20 
and the requirements of the specific application. Some of 
these different approaches are discussed below. These in no 
way represent a complete set of techniques; rather they 
highlight some key methods. The parameter solution found 
by any of these method would be well-conditioned due to the 25 
geometric diversity achieved by the LEOS, based on the 
reasoning of equ. (49). 

VI.C.l Maximum Likelihood Update 

We assume now that our estimated parameter set 0 is near 3Q 
enough the true parameter solution that there are negligible 
errors due to higher order terms, incurred in linearizing the 
measurement equation to derive (37). Based thereupon, we 
seek the maximum likelihood update: 


A 0ml = tirg max ProblAK I A©} 

QeR 4/V + S-l 1 1 J 


35 

(50) 


Consequently, the batch covariance matrix C will have 
structure: 


Cn + Crhl) Ceih , h) C e (h-> h) 
C e {t 2 ,h) C n + C e (t 2 ) C e {t 2 , r 3 ) ... 

C e (t3,h ) C e (t 3 ,t 2 ) C n + C e (t 3 ) 


(56) 


where 


V 2 ei(ti, tj ) 


C e (t h tj) = 


0 


C n 


°nl 0 

0 vis 


0 


visits tj) 


(57) 


(58) 


Given matrix C, the ML parameter update is: 




(59) 


This ML update requires knowledge of the measurement 
covariance matrix. Since the ephemeris errors are strongly 40 
correlated, the covariance matrix C=E{VV ;r } has non- 
diagonal structure. Over the course of 5 minutes of tracking, 
one can roughly model the satellite position generated from 
the ephemeris data as 

45 

r s (t)=r s (t)+Ar s (51) 

where Ar s is a constant offset, modeling the difference 
between a satellite’s true position and the position estimate 
based on the ephemeris data. We describe this offset vector 50 
in terms of normally distributed components: 

Ar=[Ax Ay Az] r ,Ax, Ay, A z ~N( 0,a e 2 ). (52) 

This error in the ephemeris data would result in an 55 

ephemeris disturbance — see equ. (40): 

e s {t) = (Psuit) - Psrit)) - ( p su {t ) - p sr (t)) (53) 

A^r) 7 ^^) + Ar^f) 7 ^) (54) 60 

fsu(t) + r sr (t) 

where the approximation is achieved with a first-order 
expansion, assuming that ||ArJ| 2 «r SM ,r 5T . Using this first- 65 
order approximation, we find the second moment of the 
ephemeris disturbance statistics: 


In essence, this iterative estimation strategy is the Gauss- 
Newton technique, where we have pre -multiplied the batch 
estimation equation (37) by the whitening matrix C~ 1/2 to 
achieve the ML update. Since the measurements in equ. (32) 
are only mildly nonlinear in r M , 2-3 iterations are sufficient 
to converge to the experimental noise floor. 

VI.C.2 Least-Squares Batch Solution via Choletsky Factor- 
ization 

The technique described above for solving the ML esti- 
mation problem requires 0(N 3 ) flops. Many techniques exist 
for reducing the computation time required by exploiting the 
sparse structure of the measurement matrix H. One such 
technique is discussed in this section. In order to preserve 
sparseness, we ignore the off-diagonal terms of the batch 
covariance matrix of equ. (56), to obtain a diagonal matrix 
C. We may pre -multiply the batch measurement equation by 
the diagonal scaling matrix ( C) -1/ 2 and then solve the 
least-squares problem. Since ( C) -1/ 2 does not change the 
block structure of H, we will not explicitly show the 
premultiplication by ( C )~ Vz , or equivalently assume the C is 
simply the identity, 1^. We may then solve the least-squares 
problem of equ. (46) by solving 

h t ha@=h t ay (60) 
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The matrix A=H 7 H has the block structure 


where 


o 

A 2,2 

0 


0 

^3,3 


A,/V+l 

^2,/V+l 

^3,/V+l 


(61) E{W(t k )W(t k ) T } = Q(t k ) 


5 


= lim 

(Ttj-KXJ 


L 4 cr^ 0 
0 0 


(65) 


M/V+1,1 An+ 1,2 ^/V+1,3 ■■■ ^/V+l,/V+l J 

where the submatrices {A, ,}, i<N+l are 4x4, {A^ +1 1 } are 
(S-l)x4 and A N+1> ^ +1 is (S-l)x(S-l). We seek a lower 
triangular matrix, L, such that LL r =A, with structure 


This process is based upon the idea that the integers are 
constant for all phase measurements and that no assumption 
is made about the user’s motion between measurements. The 
phase prediction error of equ. (36) we model with a linear 
approximation: 


L u 0 0 

0 Ei,2 o • • ■ o 


(62) AY^HiQAXiQ+Vit,) 

15 

where 


( 66 ) 


: : P/v,/v : 

.^/V+ 1,1 ^/v+i ,2 ■■■ L n+ i ; /v L/v+i,/v+i . 

20 

where the submatrices {L m }, i<N+l are 4x4 lower 
triangular, {L^ +1 ,-} are (S-l)x4 and L N+1 ^ +1 is (S-l)x(S- 
1) lower triangular. This matrix can be found via Choletsky 
Block Factorization, which is achieved with the following 2 5 
algorithm: 


o-ei(t k ) 2 + (rij 


E{V(t k )V(t k ) T } = R(t k ) = 


0 


0 1 (67) 


O-eS(tk) 2 + vis 


Denoting the covariance of our estimate AX(t^) as P(Q, 
the Kalman filtering equations for this system are: 

p(y- 1 =(p(^- 1 )+e(^ 1 ))- 1 +i/(f Jfc ) 7 ^(y- 1 ^(^)A 


for j = 1 , N + 1 


( 68 ) 


for i= j, N + 1 

l-i 
k = 1 

if i = j 

compute by Choletsky factorization Ljj s.t. LjjL J , = S 
else 

solve EijlJj j = S for L;j 


If phase -lock is attained on a new satellite s during 
parameter estimation, the initial covariance of the parameter 
estimate A“N S . is very large; this can cause computational 
difficulties. Consequently, we use instead an information 
form of the Kalman filter where we define an information 
matrix 


S(t) = P(r)- 1 = 


S Q (t) 

s qw^ T 


S N (t) 


(69) 


end 

end 

end 


and an information vector Z(t)=S(t)AX(t). Applying these 
definitions to equations (65) and (68), it is straightforward to 
show that the update equations become: 


Once L is found, A0 can be found by block back- 
substitution. 

VI.C.3 Iterative Information Smoothing 

Despite the computational efficiency, sparse matrix batch 
algorithms are difficult to administer when satellites are 
coming in/out of view, or cycle slips occur, while data is 
being stacked for processing. Information smoothing might 
be selected for its flexibility in dealing with such situations. 
In essence, the information smoother passes the linear 
Kalman filter forward and backwards over the data before 
updating the parameter set. We describe the set of param- 
eters as 


A X(t) = 


A 0(f) 
AN 


(63) 


45 


50 


S(tt) = 


Z(t k ) = 


0 0 
0 S/vA/t-i) - i S QN (t k _i) 

H(t k ) T R(t k y l H(t k ) 

0 

Zn fl/t-i ) ~ Soiyih-i) ^©(fy-i) ZQj;j-(t k -i) 

H(t k ) T R(t k r l AY(t k ) 


(70) 


In order to emulate the batch solution, the filter is passed 
55 forward and backward over the data. For the backward pass, 
we simply interchange the t^ and t^ in equ. (70), and start 
with initial conditions S(bv)=0 and Z(W)=0. The S(t^) and 
Z(t^) from forward and backward passes are linearly com- 
bined. 

60 

(71) 


where AN=[AN 2 . . . AN 5 ] r . We model the evolution of 
parameters between phase measurements as a Gauss- 
Markow process 

(64) 


Z{t k )=Z{t k ) F +Z{t k ) B (72) 

and the parameter updates are then found according to 
65 AX(t^)=S(t^) _1 Z(t^). The integer updates are found from the 
relevant vectors in AX(t a ) or AX^). This process is 
repeated until the elements of AX(t*) are negligibly small. 


NX{tk)=mk-i)+w{h) 
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VI.C.4 Solving the Max. Likelihood Problem with Integer 
Parameters 

All of the techniques discussed thus far treat the integer 
parameters as real numbers in the estimation strategy. In this 
section, we discuss a technique for solving the least squares 
problem assuming an integer parameter set. Assume that our 
current parameter estimates are close to the ML solution, so 
that no iteration is necessary and we drop the A notation. 
Consider the batch measurement of equ. (37) rewritten as 

Y=H e 6+H z z+V (73) 

where H e and are the estimation matrices for the real and 
integer parameters respectively, 0 is the matrix of real 
parameters, z is qxl matrix of integer parameters, and V is 
the disturbance matrix with covariance C. Imagine equ. (73) 
is premultiplied by the whitening matrix C -% , and the 
least-squares solution is found, with the approximation that 
z contain real elements. The resulting estimate of z can be 
shown to be normally distributed, z~N(z,<D) where 

We can consider the least-square solution z to be gener- 
ated by z=z+u where u is the disturbance term. Multiplying 
by the whitening matrix 

G = <T2, 


28 

By redefining the generation matrix G«-GF, we tighten 
the bound of equ. (76). The sum-of -squares of q independent 
normally distributed unit variance random variables is a % 2 
distribution with q degrees of freedom. We denote by 
5 F y 2(% 2 ;q) the cumulative distribution function of a % 2 ran- 
dom variable with q degrees of freedom. Once we have a 
bound on d min , d, we can find a lower bound on the 
probability of correct integer selection: 


10 


Prob{z M L = z}> F x 2 



(77) 


15 We return now to the integer least-squares estimation 
problem of equ. (75). Replacing G by G, the expression can 
be rewritten as 

Zml = org min (z - zf P^ 1 (z - Z) C8) 

20 ze2? xl 


where z=G 1 z and P^G^G) 1 . If P is diagonal, then the 
expression becomes 


Zml = arg 


min 
ze2? x 1 


I 

1=1 


(Zj - Zjf 

Pii 


(79) 


we have z=Gz+u where u~N(0, I ). The maximum likeli- 30 
hood estimate of z is then 

Zml = arg min \\z - Gz\\ 2 (75) 

zg 2? x1 


in which case we can simply find the integers by rounding: 
z MLi=\ z i\- Since G is almost orthogonal, we can use the 
rounding of z as an initial subop timal estimate of z ML , z sub . 
Since we have the lower bound 


G forms the basis, or generation matrix, of a lattice 
L(G)={Gz|zeZ <?x1 } which is conceptually illustrated in FIG. 

11 for q=2. Given the lattice L(G), 108, we can find an 
integer- valued matrix F which maps F:Z ?xl ^Z^ xl such that 
|det(F)|=l and L(GF)=L(G), 108. Given the matrix G, the 40 
algorithm of Lenstra, Lenstra and Lovasz (LLL) can be used 
to find an F such that the generation matrix G=GF has two 
key properties: 

The columns of G, {g,-}, are almost orthogonal, namely 45 

A, l 

gj = 2jWi, PM = L \Pji\ < 2 for 1 * J 


where {g,-*} are the columns of G*, the Gramm-Schmidt 
orthogonalization of G. 

|| || 2 norms of the collumns of G are bounded, namely Hg-J^ 

. . . ||g ? || 2 S2^- 1 ) /4 |det(G)|and 

Before actually determining the integer solution, we 55 
desire a lower bound on the Prob(z Aa =z). To determine this, 
we seek the radius 107, 

^min 

~~2 _ ’ 60 

of the largest ball 106 that fits within a Voroni cell 109 of 
lattice L 108. If G is orthogonalized using the Gram-Schmidt 
algorithm to G*=[g* a g* 2 . . . g* ], a lower bound on d min can 
be found: ' 65 



from equ. (76), we know that if ||z-Gz J 2 <d mi> /2 then z sub 
is the global minimum z ML . Consequently, an efficient 
algorithm for attaining the global minimum is as follows: 

* Perform the LLL on the collumns of G to generate a 

unimodular matrix F and create a new lattice generator 
matrix G=GF which is almost orthogonal. 

Zsub <- F\~G~ l z\. If II z- Gzsubh < then Zml = ^ st °p- 

* Using the algorithm of Babai [22], find a new z sub , such 

that 

q 1 
Z - GZs U b = Yj X ‘Si with | A; | < - . 

if \\z - Gzsubh < ~y~ then Zml = Zsub 311(1 st0 P- 

*tH|z-Gz„J 2 and z*^F" 1 z i „ fc 
^ while 1 

^Search for an integral point inside an ellipsoid ||z— 
Gz|| 2 <t. If no such point exists z^-e-Fz* and stop. 
*Let z* be the integral point found in the previous step. 
x«-||z . . . Gz*|| 2 . 


^ i «= m in{||g* 1 || 25 ll£* 2 ll 2 - ■ ■ \\g* q M 


(76) 
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If rs ■ 


then zml <- F£ and stop. 


1 

h~ -Pgs(ts) 


*end. 

Once the parameter {N s } are identified using any of the 
techniques discussed above, they can be regarded as con- 
stant biases. The receiver uses these estimates in equ. (35) to 
construct prediction error estimates with which to update io 
position estimates with centimeter-level precision in real 
time. One straightforward method of doing this is 

A ( 80 ) 

15 

VII. IMPLEMENTING THE INVENTION WITH 
DIFFERENT COMMUNICATION 
CONFIGURATIONS 


describes the time at which the transmission was made; and 
represents all the phase disturbances on the path from 
ground to satellite, which will almost cancel out in the 
differential measurement. The satellite downconverter mixes 
the incident signal with another at frequency oo^-oo s . The 
phase of the down-converted signal which the satellite 
transmits is then 

(j) g 

%(t s ) = (o s [t s + r s (r s )] - — p gs (t s ) - M g T s (t s ) + v gs 

The phase of the satellite signal output from the LNA 22 
of the user’s receiver 4 at time t u is then approximately 


VILA. Satellites with Multiple-Beam Downlinks 

20 

In section (VI), it was assumed that each satellite footprint 
was constructed with a single beam so that continuous 
carrier-phase data could be accumulated over the duration of 
tracking. However, rather than transmit a single beam for the 
satellite’s footprint, as in the GPS case, many LEOS trans- 2 s 
mit multiple beams. As shown in FIG. 12 which roughly 
resembles the Globalstar downlink, each beam 19a-d covers 
a portion of the satellite footprint 78. We may assume that 
each beam 79 a-d is modulated differently, and has a phase 
offset relative to adjacent beams. Imagine a receiver moves 30 
from beam a to beam b 80. There arises an interval of a few 
seconds between reaching the 3 dB-down point of beam b 81 
and the 3 dB-down point of beam a 82. During this interval, 
the receiver simultaneously tracks both beams. Before han- 
dover from beam a to beam b at time t 0 , the receiver 35 
calculates the phase different between the two beams, 
4u0o) — • The receiver then adds this difference to the 
phase measured on beam b at some later time t 1 . The sum 
then becomes which is corrected for 

any offset between the beams. 

VII.B. Compensating for the effects of a bent-pipe commu- 
nication payload 

We have assumed in section (VI) that the signal arises 
onboard the satellite. However, for a bent-pipe communica- 45 
tion payload, the satellite actually downconverts and then 
retransmits signals received from ground-based uplink sta- 
tions. This does not change the conceptual approach, but 
simply requires that additional terms be taken into account 
in the estimation strategy. An example of a bent -pipe system 50 
configuration, such as that of Globalstar, is illustrates in FIG. 

13. Consider the phase of the signal 84 a,b generated at the 
ground terminal 83 a,b: 

( 81 ) 55 


Vsutfu) = - ^p su [t u - Psu ^ )j 4 


( 84 ) 


. , / Psu(t u )\ M g ( PsuVu)\ 

{U> s - OJg)T s [t u - - —Pgs[t u - — - ) 


Psu(tu)\ 


where v su contains all the phase disturbances on the satellite - 
to-user path. Consequently, we find the user’s phase mea- 
surement corresponding to equ. (16). 


i O s ( Psu(t u )\ 

<psu{tu) = M s T u (t u ) + — Psultu - H 

c v c / 

Psu(tu)\ Mg Psu(tu)\ 


( 85 ) 


(u g - w s )r s (r„ - 
- v 5 „ - N su 2 tt 


M , Mg ( p su {tu)\ 

■) tM " ~ J - 


Similarly, the phase measurement made at the reference at 
time 4 is 


M s ( Psr(tr ) \ 

Mr) = M s T r (t r ) + —pJtr + 

c v c ) 

(Mg - W s )r s (r r - + ^-Pgs^r 


- N sr 2n 


Psr(t r ) \ 
C ! 


(86) 


As discussed in section (VI .A), the user matches the tags 
on the measurements, and performs the single difference: 

m=4su(Q-Mr) (87) 

Using a Taylor expansion and discarding all insignificant 
high-order terms, as in section (VI.A), we can recast the 
measurement in terms of the assigned timetags. The result- 
ant representation of the measurement, corresponding to 
equ. (20) is: 


We will ignore x g since its effect on the differential phase 
measurement is negligible. Using similar reasoning to that 
of section (VI.A), the phase of the uplink signal 84 a,b 
incident at the satellite 2a, b at time t s can be described: 60 

( 1 ) ( 82 ) 

= k — ~Pgs(ts) I + Vgs 


65 

Where p^ s (t s ) is the distance from the ground terminal to 
the satellite at time t s ; the expression 


, ,, m s / p su (t)^ 
Ms dp su (t ) 


s ( Psr(t) \ 

-M r -— )- 


( 88 ) 


c dt 
M s [r u (t) - r r (r)] 


[ At„ (r) + T u ] + ^ [ A 7 r (r) + T r ] 


dt 


M g dp gs (t) 


-[Ar u (t) + T u -AT r (t)-T r ]- 


c dt 
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-continued 

0J g d p gs (t) I p su (t) _ p sr (t ) 

c dt v c c 

N sr 2 jt 


- N su 2n 4 


h s (t) = 



i dpjt) 

c dt 


u g dp g s(t) 
oj s c dt 


(93) 


We will consider those terms in the expression which have 
been introduced by the bent-pipe architecture and cause equ. 
(88) to differ from equ. (20). Consider the terms 


, J / Psu(t)\ ( Psr(t)x\ 

K - u> s )[r s [t - —— J - r s [t - — jj 


(89) 


involving the instabilities of the satellite oscillators. All that 15 
has changed from equ. (20) is the multiplication factor, 
which is now (co^-co s ) instead of -co s . Therefore, the treat- 
ment of this term is similar to that of equ. (23). 

Consider the terms in equ. (88) involving 20 


dPgs 
dt ' 


We rewrite these terms as: 


25 


dPgs(t ) 

c dt 


[r«(0-r,(0] + 


a^dpgsit) 
c dt 


[r u ~ r r ] ~ 


(90) 


dp gs (t) [r T | u g dp gs (t) , p su (t) 30 

c dt u r c dt v c c / 


The 2 nd can be directly calculated if the position of the 
ground station and the satellite ephemerides are known; 35 
hence it can be actively subtracted out of the measurement. 
The 3 r d term can be ignored, using similar argruments to 
those of section (VI. A). The 4 th term is directly calculable if 
the receiver knows the satellite ephemeris, the location of 
the ground uplink, and the rough position of the user using 40 
code -phase measurements. Consequently, the fourth term is 
also calculated and actively subtracted out of the differential 
measurement. Only the 1 st term must be directly estimated. 

Eliminating all terms from equ. (88) which are either 45 
negligible, or actively subtracted out, the estimate of the 
resulting measurement corresponding to equ. (24) is 


Similarly, the estimate of our prediction error, corre- 
sponding to equ. (35) becomes 


A*M = p su [t 


(94) 


( 1 

1 - - 

V c 


dKU) dp ‘ M -)c% 0-frA-WD 


dt 


(x) s c dt 


and we may proceed with the estimation as described in 
section (VI.B). 

VII. C. Unstable Oscillators: Calibrating the LEO Oscillator 
using Navstar Satellites 

In this section, we describe how the frequency offset, or 
clock offset rate of the LEO oscillator can be calibrated 
using a GPS signal. This algorithm is only necessary for 
oscillators that have long-term frequency stabilities of 1:10 6 
or worse. Some of the mathematical steps are closely related 
to those described in detail above, so these stages have been 
left out of the explanation. The technique described here has 
been designed to be implemented completely with software, 
and requires no additional front-end hardware in the 
receiver. We will assume a bent-pipe communication archi- 
tecture for generality; the additional terms can simply be 
dropped for simpler systems. 

We can describe the phase measurements made for a 
bent-pipe LEO satellite L, 2a-c, and a Navstar satellite N 
la-d at the reference receiver 3 as: 

(95) 

(96) 

In order to cancel out the error due to the receiver’s 
oscillator 66 drift, the microprocessor 56 performs a 
weighted difference between the phase of the two satellite 
signals to find a calibration phase 


0chr) — 0Lr(^r) " 


~<pNr(t r) 


(97) 


/ Psu(t)\ Us ( Psr(t)\ 

<P s (t) * —Psu{t-—)-—Psr[t- — ) + 

1 dp ST (t) (x) g dp gs (t ) 


(91) 


50 


u s \l- 


dt (x) s c d 

{N su - N sr )2n + (p su - p sr ) + (n s 


^ j[r H (r) - T>(r)] - 

ttsr) 


55 

Following the steps outlined in section (VI .A) and section 
(VI.B), we find the measurement approximation correspond- 
ing to equ. (32): 


The incident phase from each of the the satellites may be 
described: 


, , . ( 1 / PLt(?t)\ ) U g / PLr(h)\ 

“ “ PLt (*Y " JJ “ ~ PgL^r " ) 3 

(0) L - OJg)T L {t T - PLT ^ j 


(98) 


- Plt ~ n L T 


i, , \ ( 1 / PNr(tr) / PNr(h)\ 

^ Nri^r) — — ~PNryr ~ JJ 3 " “ J “ 


PNt ~ n N t 


( Psu(t)\ ( Psr(t)\ 

y s (t) = p su {t ~ — ) - Psr[t - — ) + 


( 1 dp su (t) 

{ C dt 


Qjg dp gs (t ) 

o) s c dt 


cr(r) - N S X S 


+ n s 


(92) 60 

where the subscript g refers to the ground uplink station 
83 a,b, discussed in section (VII.B). The resultant expression 
for the weighted difference is found to be: 


65 / PNrih)\ U L / PNrih)\ (99) 

The observation matrix for the time-varying parameters, ~ (x)lTn v ~ c ) ~ ~ PNt V t ~ c ) + 

corresponding to equ. (34), becomes 
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-continued 

( PLt (*r) \ M g ( PLt hr) \ 
- p Lr{tr — ) + ~P 8 L[t T - + 

{oj 8 - w L ) r L ( f r - — — ) - N Lt 2ji + — N Nt 2jt 
\ c i ojn 


m l m l 

PLt PNt + n Lr n Nr 

(x) N 0J N 


34 

guarantee that cycle slips do not occur between bursts. If A 
is the Allen variance of limiting oscillator, we require 

Au s T s (102) 

— — <<c 1 


For Ka band downlinks for which 


The 2 nd , 2 rd and 4 th terms are directly calculated and io 
subtracted from the measurement by the reference receiver 
which knows the location of the ground uplink station, as 
well as that of the LEO and Navstar satellite. After subtract- 
ing out these terms, the calibration measurement becomes: 


<p c {t r ) = {OJ g - (x) L )T L \t T - PLT ^ j + U) L r N {t r - 


PNt hr) \ 


(jl>l (jl>l 

N Lt 2tt + — N Nt 2tt + (i Lr (i Nr 4- n Lr - 

U) N (x) N 


m l 


15 

( 100 ) 


20 


For the purpose of this calibration, we consider the LEO’s 
clock offset, T l , to be a linear function of time as a result of 
a frequency offset in the satellite oscillator. Since the GPS 25 
clock is atomic, we regard as constant. The reference 
receiver calculates the change in calibration phase A$ c over 
an interval of roughly one second, At, to find 


A t 


i (OJ g - OJ L ) 


dr L 

dt 


(101) 30 


from which 


dr L 

dt 

can be calculated with good accuracy. 

VII.D. LEO Satellites using TDMA Downlinks 40 

It should be noted that the fundamental technique of 
augmenting GPS with LEOS for geometric diversity is 
equally applicable to TDMA downlinks, where the LEO 
satellite signals 6a, b arrive in bursts of a few //sec. The time 
from the start of one burst to the start of the next is termed 45 
the scan period, T s . The time during of each burst is termed 
the receive time, T r . It has been demonstrated that continu- 
ous carrier-phase tracking of GPS C/A code-type signals can 
be achieved with high integrity for T r =2 msec, and T s =12 
msec. Although LEOS signals show more Doppler shift than 50 
GPS signals, it is well known in the art that a 3 rd order 
phase-locked loop can be implemented to maintain a run- 
ning estimate of phase (|), phase rate 



and phase acceleration 

d 2 <f> 

dF" 


* 30 GHz, 

2 n 

and T s ^25 msec, we require A «1.3xl0 -9 . The second 
limitation involves aliasing due to the dynamics of the 
receiver platform. For T s 32 25 msec, for example, the 
highest frequency component of the platform dynamics 
should not exceed half of the corresponding sampling 
frequency, or 20 Hz. Neither of these constraints is restric- 
tive for the majority of applications. 

VILE. Non-GPS Navigation Signals 

In our discussion of the carrier-phase positioning 
algorithm, we have assumed that LEOS are used to augment 
the Navstar GPS satellite fleet. It should be noted that 
Navstar GPS satellites have utility to the extend that they 
provide: 

Additional carrier-phase signal sources, which in essence 
provide additional equations so that the parameter- 
estimation problem is over-determined. 

Code-phase navigation signals which allow for correla- 
tion of the user and reference receiver clocks to better 
than ~1 //sec and also enable the receivers to achieve 
initial position estimates accurate to the meter-level. 
Other navigation satellites exist, and are planned, which 
could also fulfill both of these functions. Such systems 
include Russia’s Glonass, and Europe’s proposed Satellite 
Civilian Navigation System which will include a LEO 
segment with highly stable downlink frequencies. Our sys- 
tem of exploiting LEO satellites for geometric diversity is 
equally applicable to these other navigation satellites. For 
the general case of a mobile user, 4 such navigation signals 
should be available. If another means is used to initially 
synchronize the reference and receiver clocks, specialized 
code phase navigation signals become unnecessary. 

The technique described in section(VI) can resolve inte- 
ger cycle ambiguities on carrier-phase signals for a mobile 
user so long as a total of 5 satellites are in view. Ideally, to 
rapidly constrain all the degrees of freedom in the position- 
ing problem, 2 or more of these satellites should be in Low 
Earth Orbit. 

VII .F. Non-Differential Position Estimates 

The essential technique of augmenting GPS with LEOS 
for geometric diversity is equally applicable to the non- 
differential setting. To examine how a user receiver might 
proceed with non-differential position estimation, consider 
the first-order expansion of the phase measurement in equ. 
(18) 


OJ s ( p su (t)\ 

<psu{tu) * — pjt - + UsTuit) - 

C \ C / 

) - N su 2jt + ji su 


Ms dp suit) 
c dt 


lAW + T u ] - 


Hence, the change in phase due to satellite motion can be 
estimated over T s to maintain phase lock. Two fundamental 
limitations on the technique exist. The first concerns the 
stability of the satellite and receiver oscillators required to 


65 Since there is no need to initially estimate and compensate 
for clock offset, T M has been dropped from the expression. 
The term 
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<*>s dp su (t) T 

c dt u 


can be ignored with distance -equivalent error of <1 mm for 5 
|Tj<0.1 /usee and satellites at 1400 km. Consider now 

/ Psuit ) 

TsU 

V c 

10 

For highly stable satellite oscillators, such as atomic 
oscillators, the user can be conveyed the necessary infor- 
mation to model the clock term. Alternatively, a highly 
stable clock can be modeled linearly as T s (t)=T so +T fl t. 
Using the clock-calibration technique described in section 
(VII. C), the user may directly estimate and subtract out the 
term 


Using Ionospheric/Tropospheric models and/or dual fre- 
quency ionospheric calibration, the user receiver could esti- 
mate and largely subtract out the error terms n su , leaving a 25 
residual An su . Converting the resultant estimate to distance - 
equivalent form, 


, x / Psu(tu)\ {, 1 dPsuitu)) , . 

) + { l ~c ~~d^r u{tu) " 


CTso - X S N SU + + kn su 


(103) 

30 


Using one of the Navstar satellites, say satellite 1, as a 
reference satellite for differencing, we make similar redefi- 35 
nitions to those described in section (VI.B) 


y su (t u ) = yM + A hu 

at 




Tji 


r u (t u ) = r a (rj ■ 


Ai 


X s ) \ s \ 


N lu -N u 


2n ) 


Pi u \ 
' 2^1 


(104) 


40 


For the _Navstar satellites, for which T, o =0, the new 45 
parameter N su remains integer- valued. The measurement 
may then be described 




Psu[tu - 


Psu(tu)\ , (, 1 dPsuitu)) _ . . 


(105) 

50 


X S N SU + A n su 


One skilled in the art could then proceed with a similar 
estimation strategy to that described in section (VI.B), with 55 
the positioning accuracy depending primarily on the mag- 
nitude of An SM . 

VIII. EXAMPLE IMPLEMENTATION OF THE 
INVENTION: AUGMENTING GPS WITH THE 
GLOBALSTAR TELECOMMUNICATIONS 60 

CONSTELLATION 

VIII. A. The necessary criteria for cm-level positioning 
In order for a LEO constellation to rapidly resolve cycle 
ambiguities with integrity for a mobile user, the following 
criteria should be fulfilled: 65 

There should ideally be 2 or more LEO satellites available 
for tracking. 


36 

A carrier signal should be traceable for a time period of a 
few minutes. 

The satellite ephemerides should be known to good 
accuracy. 

The SNR ratios should be sufficient for accurate carrier 
phase estimation. 

All of these criteria are fulfilled by the Globalstar Con- 
stellation. Carrier phase from one satellite can be tracked for 
several minutes at a time. In addition, the constellation has 
GPS sensors onboard so the ephemerides can be estimated 
to <20 m rms. FIG. 14 describes the percentage availability 
of the satellites at different latitudes. Note that there are 
always 2 satellites available above 10*** elevation over the 
continental United States. 

VIII. B. Integrity with RAIM 

In all availability and performance analyses, we assume 
that the LEOS are functional, and that no cycle slips occur 
over the tracking duration. For high-integrity applications, 
this assumption cannot be made, and position solutions 
would be independently validated via receiver autonomous 
integrity monitoring (RAIM). In essence, the RAIM algo- 
rithm checks if the residual of the least-squares solution at 
each time t, ||AY(t)-H(t)A0(t)||, is greater than some thresh- 
old R. R is set to meet a continuity requirement — that is, not 
to exceed the allowed number of false alarms of system 
malfunction caused by regular measurement noise. For a 
given acceptable error radius a, we can only guarantee that 
RAIM will alert us to position errors ||? M (t)-r M (t)||>a using 
the threshold R for given satellite geometries. FIG. 15 
addresses the availability of such geometries with respect to 
latitude to alert for radial errors of 1.1 m while allowing a 
continuity risk of 2xl0 -6 per 15 sec. These results were for 
122.17** West. We assume that GPS is augmented only 
with the Globalstar Constellation. A conservative phase 
noise variance of 1.4 cm, and phase reading rate of 5 Hz are 
assumed. 

VIII. C. A Joint GPS/Globalstar Reference Transceiver 

Each Globalstar satellite downlink contains 13 1.23 MHz 
bands, spanning 2483.5 MHz to 2500 MHz in the S-band. 
Each band supports 128 CDMA channels, one of which 
harbors a pilot signal, the modulation of which can be 
described by equ. (1). In this context, D(t) refers to an outer 
PN sequence of length 288 which is chipped at 1.2K cps. 
C,(t) is created from the sum of 2 inner PN sequences, each 
of length 2 10 , which is then filtered to a 1.23 MHz band- 
width. C e (t) is a different code, created in a similar manner 
to C,(t). FIG. 16 illustrates an architecture for the joint GPS 
and Globalstar reference receiver 3, together with a trans- 
mitter subsection 90. The rf filter 23 has a center frequency 
of f G =2492 MHz and a bandwidth of roughly 75 MHz. The 
^Grf 24 and f Gif 28 mixing frequencies can be generated 
using an integer-N synthesizer 89 with a dual modulus 
prescaler. The if mixing scheme 27 adheres to that of FIG. 
4b. The signal is filtered 36 to a pre -correlation two-sided 
bandwidth, B c of less than «2.5 MHz. In-phase and quadra- 
ture sampling 37 of the Globalstar signal occurs at f s «20 
MHz, the clock speed 99 of the Globalstar tracking assembly 
88 circuitry. In the preferred embodiment, the LRU 8 is 
implemented using a ground-based transmitter 90 which 
employs a VHF rf carrier frequency f r The f f carrier signal 
96 is generated by the frequency synthesizer 89. The data 
modulation section 91 is controlled by the microprocessor 
56 via the data/control bus 93. Data is modulated onto the 
carrier using PSK modulation at a rate of 9600 baud. The 
signal is then amplified 94 and transmitted via a VHF 
antenna 95. 
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VIII. D. Primary Error Sources for Precise Navigation 
This section gauges those additional sources of error 
which cannot be directly estimated by the data reduction 
algorithm of section VI.A. We provide only enough detail to 
roughly calibrate the sources of error. 

VIII.D.l Receiver Phase -Tracking Errors 

We will gauge this error using equ. (9), since the filtering 
of the Globalstar inner PN sequences does not have a 
substantial effect on the phase tracking performance. It is 
prudent, for most implementations of the invention, to select 
a narrow phase-locked loop bandwidth B L ~10 Hz. This B L 
enables a second-order phase-locked loop with damping 
ratio £«0.7 to track phase acceleration of 100.2jtrad/sec 2 
with an error <0.1 rad. For a receiver noise figure of roughly 
3 dB, the nominal Globalstar transmission archives 

A 2 

- — = 37.5 dB - Hz. At B L = 10 Hz, 


we expect a 1-a phase error due to thermal noise or roughly 20 

0.12 rad. 

VI II. D. 2 Ionospheric Errors 

The distance -equivalent group delay due to the iono- 
sphere can be as large 20 m. However, if a differential 
carrier-phase measurement is taken as in equ. (20), and the 25 
user and reference are within 10 km, the resultant error is 
governed by local irregularities in the ionospheric structure, 
which delay the signal to the user and reference station by 
different amounts. For S-band transmissions, and a user- 
reference separation of d km, we estimate the resultant phase 30 
error as a normally distributed random process of zero mean 
and variance X 2 4.4xl0 _4 d rad 2 , where X is the wavelength in 
cm. This leads to 1-a phase deviations on Globalstar signals 
of 0.25 rad and 0.57 rad for distances of 1 km and 5 km 
respectively. The corresponding GPS deviations are 0.44 rad 35 
and 0.99 rad respectively. 

VI II. D. 3 Tropospheric Errors 

Without any form of differential correction, the delays 
cause by the troposphere at a satellite elevation of 10t deg are 
roughly 14 m. With differential measurement, and a baseline 40 
separation of d km<10 km the remaining tropospheric delay 
can be roughly modeled as a normal distribution ~N(0 
cm, (0.1 d) 2 cm 2 ). 

VIII. D. 4 Ephemeris Errors 

For the Globalstar satellites, we expect the ephemeris 45 
disturbances discussed in section (VI.A) to be bounded by 
1.5 cm and 7.2 cm for d of 1 km and 5 km respectively. The 
ephemerides of the GPS satellites are known to within 
roughly 10 m rms; the resulting ephemeris errors are 
bounded by 0.05 cm and 0.25 cm respectively. 50 

VIII. E. Expected performance of a System Using only the 
Globalstar Constellation 

To illustrate the performance of the invention, we discuss 
a Monte-Carlo simulation which indicates the behavior of a 
system augmenting GPS with the Globalstar Satellites. 55 
Separate simulations are conducted for separations between 
user and reference station of 1 km and 5 km. The simulations 
are all conducted assuming the user is in Palo Alto, Calif., 
and is capable of seeing all satellites above 10*** elevation. 

For this study, conservative disturbances were assumed due 60 
to ionospheric and tropospheric delay, thermal noise, 
ephemeris errors, and oscillator instabilities. For each 
simulation, the assumed time of the experiment was varied, 
sequentially sampling the interval of 12 hours, roughly one 
period of the Navstar satellite orbits. For each simulation, 65 
the user receiver is ascribed velocity and displacement from 
the reference receiver in a random direction. It was assumed 
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that the user’s motion was relatively slow, so that the 
user-reference separation was roughly constant over the 
course of tracking. 

FIG. 17 displays the 1-a deviation of radial position 
errors as a function of tracking time for a mobile user. The 
parameter estimates were found using equ. (59). Each point 
corresponds to the mean error deviation averaged over 200 
simulations. FIG. 18 displays the evolution of the lower 
bound on the probability of selecting the correct set of 
integers for the Navstar satellites in view over the tracking 
period. Each point in this figure represents the worst case for 
200 simulations. This lower probability bound is calculated 
according to the technique discussed in section (VI.C.4). 
Once the integers are correctly identified, a user may rely 
completely on GPS measurements were ephemeris errors are 
negligible. Since the measurement errors for each satellite 
are roughly distributed as ~N(0, a„/), the positioning error 
deviation for a static user decreases roughly as 

where N is the number of measurements taken on each 
satellite. 

“The present invention has been described in terms of a 
preferred embodiment. The invention, however, is not lim- 
ited to the embodiment depicted and described. Rather, the 
scope of the invention is defined by the appended claims.” 
What is claimed is: 

1. A method for estimating a position of a user device in 
a satellite -based navigation system, the method comprising: 

transmitting carrier signals from a set of satellites, 
wherein the set of satellites includes a set of LEO 
satellites; 

accumulating and sampling at a reference station the 
carrier signals to obtain reference carrier phase infor- 
mation comprising geometrically diverse reference car- 
rier phase information from the set of LEO satellites; 
accumulating and sampling at the user device the carrier 
signals to obtain user carrier phase information com- 
prising geometrically diverse user carrier phase infor- 
mation from the set of LEO satellites; and 
calculating the precise position of the user device based 
on the reference carrier phase information and the user 
carrier phase information, wherein the geometrically 
diverse reference carrier phase information and geo- 
metrically diverse user carrier phase information from 
the set of LEO satellites are used to resolve parameters 
related to integer cycle ambiguities in the reference 
carrier phase information and the user carrier phase 
information. 

2. The method of claim 1 further comprising: 
receiving code signals from a set of navigation satellites; 
measuring at the reference station the code signals to 

obtain reference code phase information; 
measuring at the user device the code signals to obtain 
user code phase information; 
estimating user and reference clock biases from the user 
and reference code phase information; and 
correcting for clock offsets using the estimated user and 
reference clock biases. 

3. The method of claim 1 further comprising initializing 
a device navigation algorithm by estimating an approximate 
user position using code phase signals received from a set of 
navigational satellites. 
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4. The method of claim 1 further comprising communi- 
cating differential code phase correction data and the refer- 
ence carrier phase information from the reference station to 
the user device. 

5. The method of claim 1 further comprising communi- 5 
eating LEO satellite ephemeris data to the user device 
directly from the reference station, or using a satellite data 
link. 

6. The method of claim 1 wherein the step of calculating 
the precise position of the user device comprises predicting 10 
present reference carrier phase information based on past 
reference carrier phase information. 

7. The method of claim 1 wherein the step of calculating 
the precise position of the user device comprises compen- 
sating for frequency dependent phase delay differences 15 
between navigation carrier signals and LEO carrier signals 

in user and reference receiver circuits. 

8. The method of claim 1 wherein the step of accumu- 
lating and sampling the carrier signals at the user device 
comprises reading navigation carrier information and LEO 20 
carrier information within a predetermined time interval 
selected in dependence upon an expected motion of the user 
device and the motion of the set of LEO satellites. 

9. The method of claim 1 wherein the calculating step 
comprises accounting for a carrier phase offset between two 25 
LEO beams from a single LEO satellite. 

10. The method of claim 1 wherein the calculating step 
comprises calibrating LEO oscillator instabilities using 
navigation satellite information. 

11. The method of claim 1 wherein the calculating step 30 
comprises compensating for phase disturbances resulting 
from a bent pipe LEO communication architecture. 

12. The method of claim 1 further comprising the step of 
monitoring the integrity of the calculating step. 

13. A satellite-based navigation system comprising: 35 

a set of satellites to transmit carrier signals, wherein the 

set of satellites includes a set of LEO satellites; 

a reference station to track the carrier signals to obtain 
reference carrier phase information comprising geo- 
metrically diverse reference carrier phase information 40 
from the set of LEO satellites; 

a user device including a receiver to track the carrier 
signals to obtain user carrier phase information com- 
prising geometrically diverse user carrier phase infor- 45 
mation from the set of LEO satellites and a micropro- 
cessor to calculate a precise position of the user device 
based on the reference carrier phase information and 
the user carrier phase information, wherein the micro- 
processor uses the geometrically diverse reference car- 5Q 
rier phase information and geometrically diverse user 
carrier phase information from the set of LEO satellites 
to resolve parameters related to integer cycle ambigu- 
ities in the reference carrier phase information and user 
carrier phase information; and 55 

a communications link between the reference station and 
the user device. 

14. The system of claim 13 wherein the set of satellites 
further comprises a set of navigational satellites and wherein 
the communications link conveys differential code phase 60 
correction data and the reference carrier phase information 
from the reference station to the user device. 
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15. The system of claim 13 wherein the communications 
link conveys LEO satellite ephemeris data to the user device 
from the reference station. 

16. A user device for providing satellite -based navigation, 
the device comprising: 

at least one antenna to couple to carrier signals transmit- 
ted from a set of satellites wherein the set of satellites 
includes a set of LEO satellites; 
a first receiver for to track the carrier signals to accumu- 
late and sample carrier phase information comprising 
geometrically diverse user carrier phase information 
from the LEO satellites; 

a second receiver, not necessarily distinct from the first 
receiver, to obtain reference carrier phase information 
transmitted from a reference station, wherein the ref- 
erence carrier phase information comprises geometri- 
cally diverse reference carrier phase information from 
the set of LEO satellites; and 
a microprocessor to calculate a position of the user device 
based on the reference carrier phase information and 
the user carrier phase information, wherein the micro- 
processor uses the geometrically diverse reference car- 
rier phase information and geometrically diverse user 
carrier phase information from the set of LEO satellites 
to resolve parameters related to integer cycle ambigu- 
ities in the reference carrier phase information and user 
carrier phase information. 

17. The device of claim 16 wherein: 

the set of satellites further comprises navigation satellites; 
the first receiver measures navigation code signals to 
obtain user code phase information; 
the second receiver receives reference code phase infor- 
mation transmitted from the reference station; and 
the microprocessor estimates user and reference clock 
biases from the user code phase information and ref- 
erence code phase information, and uses the estimated 
clock biases to correct for clock offset errors. 

18. The device of claim 16 wherein the first receiver reads 
navigation carrier phase information and LEO carrier phase 
information at times separated by no more than a predeter- 
mined time interval which is dependent upon expected 
movement of the device and the movement of the LEO 
satellites. 

19. A user device for providing satellite -based navigation, 
the device comprising: 

at least one antenna to couple to signals transmitted from 
a set of satellites, wherein the set of satellites includes 
a set of navigation satellites and a set of LEO satellites; 
a receiver to track the signals to obtain code phase 
information and carrier phase information comprising 
geometrically diverse carrier phase information from 
the set of LEO satellites; and 
a microprocessor to calculate a position of the user device 
based on the code phase information and the carrier 
phase information, wherein the microprocessor uses the 
geometrically diverse carrier phase information from 
the set of LEO satellites to resolve parameters related 
to integer cycle ambiguities in the carrier phase infor- 
mation from the set of navigation satellites. 
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